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Chapter 1 
INTRODUCTION 


The theoretical foundations and user instructions for a FORTAN code for the 
design and analysis of composite grid-stiffened cylinders subjected to global and local 
buckling constraints, and strength constraints are presented using a discrete optimizer 
based on a genetic algorithm. An improved smeared stiffener theory is used for the 
global analysis. Local buckling of skin segments are assessed using a Rayleigh-Ritz 
method that accounts for material anisotropy. The local buckling of stiffener segments 
are also assessed. Constraints on the membrane strains in the skin and stiffener 
segments are imposed to include strength criteria in the grid-stiffened cylinder design. 
Design variables are the axial and transverse stiffener spacings, stiffener height and 
thickness, skin laminate stacking sequence, and stiffening configuration. The design 
optimization process is adapted to identify the best suited stiffening configurations 
and pattern for grid-stiffened composite cylinder with the length and radius of the 
cylinder, the design in-plane loads, and material properties as inputs. 

The theoretical foundations for the analyses involved in the buckling of grid- 
stiffened circular cylinders are discussed briefly in Chapter 2. Instructions for setting 
up input files for the FORTRAN code are given in Chapter 3. To provide flexibilities in 
performing different types of optimization, instructions are also provided in Chapter 
3 for modifying the code to adjust the number of design variables to facilitate a 
particular type of optimization. 


2 


Chapter 2 

THEORETICAL 

FOUNDATIONS 


The buckling analysis of grid-stiffened composite circular cylinders subjected to 
combined loads requires several key steps. Herein, acceptable designs are those which 
buckle globally and do not exhibit any local skin buckling or stiffener crippling, and 
the membrane strains in the skin and stiffener segments are below an acceptable level. 
The first step in the design process is to assess the global buckling response of a grid- 
stiffened shell. Once this global buckling response is determined, the second step is 
to determine the local skin buckling response for the quadrilateral and/or triangular 
skin segments between the stiffeners. The third step is to determine whether stiffener 
buckling or stiffener crippling has occurred at this global buckling load level. Finally 
the membrane strains in the skin and stiffener segment are determined. 

The theoretical foundations for the various steps are 

• Buckling of simply-supported orthotropic cylinders. 

• Improved smeared stiffener theory ([1]). 

• Buckling of panels with general parallelogram and triangular-shaped 
planform ([2, 3]). 

• Crippling of stiffener segment ([4]). 

• Load distribution between skin and stiffeners. 

• Strain analysis. 

• Optimization strategy. 


and are discussed in this chapter. 



2.1 BUCKLING OF SIMPLY-SUPPORTED ORTHOTROPIC 
CYLINDERS 
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The global buckling analysis is based on a Rayleigh- Ritz method using a first- 
order, shear-deformation theory and the improved smeared-stiffener modeling ap- 
proach discussed in [1]. The cylinder is assumed to be simply supported and hence, 
the Rayleigh-Ritz method for the global analysis assumes the following Ritz functions 
for the axial displacement («), the circumferential displacement (u), the transverse 
displacement (in), and the cross-sectional rotations 4> x and (f) y \ 
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where L and R are the length and radius of the cylindrical shell, respectively, and N 
is the number of terms in the Fourier series. The coordinate system for the cylinder 
is shown in Figure 1. The minimum potential energy principle is used with the 
Rayleigh-Ritz method, and include Sanders-Koiter shell theory ([5, 6]). 


2.2 IMPROVED SMEARED STIFFENER THEORY 

The improved smeared stiffener theory for stiffened cylinders, used here includes 
skin-stiffener interaction effects. Skin-stiffener interaction effects may lead to overes- 
timation of buckling loads especially when the stiffener spacings are not small. 

If a stiffened plate is bent while it is supported on all four edges, the neutral 
surface in the neighborhood of the stiffener will lie between the mid-plane of the 
skin and the centroid of the stiffener. It is convenient to think of this as a shift 
of the neutral surface from the centroid of the stiffener. Hence, the approximate 
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stiffness added by a stiffener to the skin stiffness will then be due to the skin-stiffener 
combination being bent about its neutral surface rather than due to the stiffener being 
bent about its own neutral surface or the skin neutral surface. The shifted location 
of the neutral surface is determined theoretically through a study of the local stress 
distribution near the skin-stiffener interface for a panel with a blade stiffener. 

The neutral surface profile of the skin-stiffener combination is developed ana- 
lytically using the minimum potential energy principle and statics conditions. The 
skin-stiffener interaction is accounted for by computing the bending and coupling 
stiffness due to the stiffener and the skin in the skin-stiffener region about a shifted 
neutral axis at the stiffener. 

A grid-stiffened cylinder may be considered to be an assembly of repetitive units 
or unit cells (see Figure 1). Any stiffener segment in the unit cell may be isolated in 
a semi-infinite skin-stiffener model as shown in Figure 2 for a diagonal stiffener. The 
approach for obtaining the neutral surface in a semi-infinite stiffened panel is given 
in Reference [1]. 

A typical profile of the neutral surface for a skin- stiffener combination is shown in 
Figure 3. The distance y* represents the distance from the centerline of the stiffener 
to the point where the neutral surface coincides with the mid-surface of the skin. The 
average of the neutral profile over the distance y * is Z*. The quantities y* and Z* 
are obtained numerically. 

The smeared stiffnesses of a stiffened panel is obtained by mathematically con- 
verting the stiffened panel to an equivalent unstiffened panel (Ref. [7]. The smeared 
stiffnesses are developed on the basis that the strain energy of the stiffened panel 
should be the same as that of the equivalent unstiffened panel. These smeared stiff- 
nesses can then be used in a Rayleigh- Ritz type analysis to solve for buckling loads 
of the stiffened panel. In Reference [7], the strain energy of the skin and stiffeners 
in the unit cell is obtained by using stiffnesses of the skin and the stiffeners which 
are computed about the mid-surface of the skin. Since, there is a shift in the neutral 
surface at the stiffener, the stiffness of the stiffeners and the skin segment directly 
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above it has to be computed about a shifted neutral surface so as to account for the 
skin-stiffener interactions. 

The correction to the smeared stiffnesses due to the skin-stiffener interaction is 
herein introduced by computing the stiffness of the stiffener and the skin segment 
directly contiguous to it according to the following criteria. 

1. If y* < t /4 j then the reference surface for the stiffener is Z n . 

2. If y* > t/ 4, then the reference surface for the stiffener is Z* . 

In either case, the reference surface of the skin is taken to be its mid-surface. Other 
more elaborate and accurate schemes can be used to introduce the skin-stiffener inter- 
action using the neutral surface profile. However, the one described herein is simple, 
and provides sufficiently accurate buckling loads for the preliminary structural design 
([!])• 


2.3 LOCAL BUCKLING OF SKIN SEGMENTS 

The shape of a skin segment on a grid-stiffened panel depends on the stiffening 
configuration. If the stiffening configuration involves diagonal stiffener only, then the 
skin segment has a rhombic planform. If the stiffening configuration has diagonal 
stiffeners with axial or transverse stiffeners, then the skin segment has an isosceles 
triangular planform. For a general grid- stiffened panel, the skin segment has a right- 
angle triangular planform, and for an isogrid panel the skin segment has an equilateral 
triangular planform. 

Buckling analyses for panel with these kinds of planforms is achieved through the 
use of ” circulation function” and accounts for material anisotropy, different boundary 
conditions, and combined in-plane loading. A First-Order Shear- Deformation Theory 
is used. The shell theory that is used can be either Sanders-Koiter, Love, or Donnell 
theory. This is achieved through tracer coefficients. 
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2.3.1 Physical and Computational Domain 

The buckling analysis of these local skin segments is enhanced by mapping their 
physical domain into a computational domain. Consider a general quadrilateral or 
triangular panel subjected to a state of combined in-plane loading where the loading 
and material properties are defined using the coordinate system (x — y) shown Figure 
4. The transformation from a physical domain to computational domain is necessary 
when dealing with general quadrilateral and triangular geometries in order to facilitate 
the computation of linear stiffness and geometric stiffness matrices and imposition of 
boundary conditions. 

The physical domain V[x,y] is transformed to a computational domain V[^y] 
as indicated in Figure 4 The mapping for a quadrilateral is 

x (t>v) = 

i=l 

y(^v) = Y, N i(Z,y)yi ( 2 ) 

i = 1 

where X{(i = 1,2, 3,4) and yi(i = 1,2, 3, 4) are the physical coordinates of the i th cor- 
ner of the panel, £ and y are the natural coordinates for the quadrilateral geometries, 
and N{ ( i = 1,2, 3, 4) are the bilinear mapping functions given by 

N i (£>»?) = j(l -£)(! + >?) 

m,v) = J(i +e)(i+'?) 

m,v) = l(i+«(i-7) 

N 4 (Z,V) = 1(1- {)(! -ij) 


The Jacobian of the transformation is 


• dx 

di£ ' 

di 


dx 

dy 

. dri 

dr] . 


(3) 


which is independent of the natural coordinates for general parallelogram-shaped ge- 
ometries. This results in substantial computational savings in the overall formulation. 
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The mapping for a general triangle is 

z(f ? V, P ) = + PX2 + px 3 

ytf,rhP) = &/i + % + P2/3 (4) 

where £, 77 and p are the area coordinates for the case of triangular geometries, and 
Xi(i = 1,2,3) and y t (i = 1,2,3) are the physical coordinates of the i th corner of the 
panel. Note that the third area coordinate will be expressed in terms of the other 
two or p = (1 — £ — 77) based on the constraint that the sum of the area coordinates 
must be equal to one. The Jacobian of the transformation is independent of the area 
coordinates. The Jacobian, in either case, is used to relate derivatives in the two 
domains. 

2.3.2 The Rayleigh-Ritz Method 

The Rayleigh-Ritz method is an approximate method for solving a certain class of 
problems. Accordingly, trial functions with some unknown coefficients and satisfying 
the essential or geometric boundary conditions are introduced in the energy functional 
of the problem. The minimum conditions of this functional are then imposed, and 
resulting algebraic equations are solved for the unknown coefficients. These trial 
functions are called the “Ritz” functions. 

The Ritz functions used here are expressed in terms of natural coordinates for 
the quadrilateral geometry or area coordinates for the triangular geometry for dis- 
placement field. The components of the displacement vector are three translations 
(Di, D 2 , D% — uo,Vo,w) and two cross-sectional or bending rotations (D 4 , D 5 — 
when considering transverse-shear deformation effects. Each displacement 
component is approximated independently by a different Ritz function. The approx- 
imation for the i th component of the displacement vector is given by 

N 

(£ , 7 ) = 

3=1 
N 

= f°r * = 1,2,3, 4,5 

3 = 1 


( 5 ) 
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where d{j represents the j th term in the TV-term approximation for the i th displace- 
ment component, are unknown coefficients to be determined, and P t (^,7/) are the 
circulation functions. 

The circulation functions I\- in Equation (5) are then used to impose different 
boundary conditions along each edge of the plate. Each term I\ is the product of 
three functions in the case of the triangular plate geometry and four functions in the 
case of the quadrilateral plate geometry. Each function is the equation of an edge of 
the triangular or quadrilateral plate as shown in Figure 4 raised to an independent 
exponent for each displacement component. Thus, the circulation functions for the 
quadrilateral plate are 

r. = (i->?r(i-^‘(i+>?) ri (i+er 

and for the triangular plate are 

r,- = ^(i (6) 

For example, considering the quadrilateral plate case, pi refers to edge 1, qi refers 
to edge 2, r 4 - refers to edge 3, s; refers to edge 4 as indicated in Figure 4. These 
exponents are used to impose different boundary conditions. If the i th displacement 
component is free on a given edge, then the exponent for that edge will have a value 
of zero. If the i th displacement component is constrained on a given edge, then the 
exponent for that edge will have a value of one. Only geometric boundary conditions 
are imposed in this approach. Thus, a simply supported condition for bending fields 
can be imposed on edge 1 by setting: 

• pz = 1 for w, p^ = 0 for <j> x , p 5 = 0 for (f> y 

A clamped condition for bending fields can be imposed on edge 1 by setting: 

• pz = 1 for w, p 4 = 1 for <j) x , p$ — 1 for (f> y 

A free-edge condition can be imposed on edge 1 by setting: 

• pi — 0 for ^o, ^o, w , <f> x and <f) y 
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The term fj in Equation (5) is a polynomial function in £ and 7 ?, and in its 
simplest form is a power series in £ and 77 and is expressed as 

MZ>v) = rV' 

mi, n, = (0, 0), (1,0), (0, 1), (2, 0), (1, 1), (0, 2), ... (7) 

The values of rrij and rij are used basically to define terms in a two-dimensional Pas- 
cal’s triangle. The number of terms N in Equation (5) defines the order of a complete 
function in two variables. The table below gives the value of N for polynomials of 
different degrees. 

Table 1 Degree of polynomials with value of N 


N 

Degree of 
polynomials 

N 

Degree of 
polynomials 

1 

0 

45 

8 

3 

1 

55 

9 

6 

2 

66 

10 

10 

3 

78 

11 

15 

4 

91 

12 

21 

5 

105 

13 

28 

6 

120 

14 

36 

7 

136 

15 


2.4 CRIPPLING OF STIFFENER SEGMENT 


The local stiffener segment is analyzed to determine whether stiffener crippling 
will occur. Reference [4] provides a method for determining the buckling load of 
a stiffener segment. Accordingly, the stiffener segment at the nodes or intersection 
points of stiffeners are assumed to be clamped while the stiffener-skin attachment is 
assumed to be a simple support. From Ref. [4], the crippling load of the stiffener is 
Ncrip and is given by 


N, 


crip 



N cl 


3 r 47 t 2 E u Gu 

a h2Ll[l-(ui 2 E22/En)] h> 


where 


( 8 ) 
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where 

s z = |Gi3^ s , is a shear correction factor, 

L\ = 2 L is the length of the stiffener, 
h is the height of the stiffener, 
h is the height of the stiffener, 
t s is the thickness of the stiffener, 

En is the longitudinal modulus of the stiffener material, 

E22 is the transverse modulus of the stiffener material. 

2.5 LOAD DISTRIBUTION BETWEEN SKIN AND STIFFENERS 

The global buckling load is assumed to be a scalar multiple of the design load 
and has the form 

(N x , N y N xy ) = Xg (Ah, N 2 N u ) (9) 

where Ah, Ah, and Ah2 are the applied in-plane prebuckling loads and represents the 
design load. Once the global buckling load factor (Ag) has been determined using the 
improved smeared stiffener theory, the loads acting on the stiffener and skin segments 
have to be determined by distributing the loads based on the extensional stiffness of 
the skin and the stiffener. 

The loads acting on the skin and stiffener segments are computed based on a 
global load factor of A q and these loads are used to determine the local buckling 
load factor of the skin, ( \ s k ), local crippling factors of axial stiffener segment, (A^, 
transverse stiffener segment, (A2) and diagonal stiffener segment, (A3). These local 
buckling and crippling load factors describe the buckling characteristics of the stiff- 
ened cylinder and is as follows 

• For A S k, Ai, A2, A 3 > 1.0, then the cylinder buckles globally at an axial load of 

Ag Ah, be., X cr — Ag- 

• If one of A s jt, Ai, A 2 , A 3 < 1.0, then the stiffened cylinder buckles locally. If \ s k 
< 1.0, then skin buckling occurs, and if Ai < 1.0 then crippling of the axial 
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stiffener occurs. For this case, A cr = A; x \q where A t - is < 1.0, and subscript 
i can be any one of sk , 1, 2 or 3. 

• If more than one of A 5 j > vy Ai, A2 and A3 are < 1.0, then local buckling of the 
stiffened cylinder occurs and A cr = A ,■ x Ag where A t - is the minimum of any of 
A sk, Ai, A 2 or A 3 with values < 1.0. 


The procedure ([8, 9] for distributing the applied loads for a general grid-stiffened 
circular cylinder are computed as follows: 

_ 2(An)i h 2(An) 3 h sin 3 6 
{Au)t = ^ 1 £ r f/iiijs 

2(y4h ) 2 h 2(An) 3 h cos 3 9 

[A 22)t — 1 r 1^22 )s 

a a 

2(An) 3 h cosO sin 2 0 

(A 66 )t = + (^6e) 5 


(10) 


where is total smeared axial extensional stiffness of the grid-stiffened panel, 

(A 22 )t is the total smeared transverse extensional stiffness of the grid-stiffened panel, 
(A 6 6 )t is the total smeared in-plane shear stiffness of the grid stiffened panel, (An)i, 
(An) 2 , (Ah) 3 are the extensional stiffness of the axial, transverse and diagonal stiff- 
eners, respectively, (A{j) s is the extensional stiffness of the skin, 0 is the orientation 
of the diagonal stiffener, and h is the height of the stiffener. Second, the loads carried 
by the skin segment which could be either a general parallelogram- shaped geometry 
or a general triangular-shaped geometry, at the panel global buckling load are 


( N x ) s k 

{Ny) S k 




(^ll)s 

(An)r 

(^22)3 

(A 2 2)r 

(^66) s 

{Aee)T 


N x = 

Ny = 
N xy = 


(^11)3 

(^-22)3 

(^22)t 

(^66)3 

(A.66)t 


A G JV, 


Ag^2 


A G N Vt 


( 11 ) 


These values then correspond to the design loads used for the in-plane prebuckling 
load in the skin-segment local buckling computation. If the critical buckling load 
factor of the skin segment \ s k is greater than or equal to one, then the skin-segment 
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buckling load is greater than or equal to the global buckling load of the grid-stiffened 
panel. Third, the loads carried by each stiffener are computed. The load carried by 
the axial stiffener is 

Wi = (12) 

where N cr i P is determined using Equation (8). The critical buckling load factor, Ai , 
of the axial stiffener has to be greater than or equal to one. The load carried by the 
transverse stiffener is 

(N x )t = 14 ^ *gN 2 = A 2 N crip (13) 

\A-22)T 

and the critical buckling load factor, A 2 , of the transverse stiffener has to be greater 
than or equal to one. The load in the diagonal stiffeners has components from the 
axial, transverse, and in-plane shear loadings and is given by 


( Nx)z — NdxSind T NdyCosO T (-/V^y^cos^ T {^Ndxy^jySinO — A3 Ncrip 


where 


Nd x = 

{A f/i n3e x aNl 

(Ah)t 

II 

£ 

oN 2 

(A 22 )T 

( N dxy ) x — 

( An) 3 cos6sin 2 9 b 

t a \ 12 

(A 66 )t a 

{Nd Xy )y — 

(A n ) 3 cosOsin 2 e 

( A x agN 12 

t^66jT 


(14) 


Ndx is the contribution from the axial in-plane loading, Nd y is the contribution from 
the transverse in-plane loading, (Nd xy ) x is the contribution from the in-plane shear 
loading along the edge where x is constant, and (Nd xy ) y is the contribution from the 
in-plane shear loading along the edge where y is constant. The critical buckling load 
factor, A 3 , of the diagonal stiffener has to be greater than or equal to one. 


2.6 STRAIN ANALYSIS 

The critical buckling load factor of the stiffened cylinder is A cr where A cr takes 
on values as discussed in Section 2.5 and based on this load value the loads acting on 
the skin and stiffeners segments are obtained. For an axial load in the skin segment 
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of Nxski a circumferential load of N ya kt a shear load of N xys k and the loads in axial, 
transverse and diagonal stiffener segments of N x i, N x 2 and N x 3 , respectively, the 
membrane strains in the skin and stiffener segments are 


,0 

xsk 

= a\\ k) N x . k 

+ 

a 

a 12 

Nysk 

+ 

J Sk ) 
a 16 

N xys k 

N 

ysk 

= a[?N x . k 

+ 

( sk ) 
«22 

Nysk 

+ 

J sk ) 

a 26 

N X ysk 

1 xysk 

~ a[t ] N x , k 

+ 

a ^ 
a 26 

Ny sk 

+ 

( sk ) 

a 66 

N X ysk 

4 

= a { "N xl 







e x2 

= atfN x2 







6 x3 

= <$N x3 








where e° sA ., e° sA; , and / y xysk - ) are the axial, circumferential, and shear membrane strains 
in the skin. e° 2 , and e° 3 are the membrane strain in the axial, transverse, and 
diagonal stiffener segments respectively. The quantities a\j k \ a^, and are 
axial flexibilities of the skin, axial, transverse, and diagonal stiffeners. 


The strain level factors for the skin, axial, transverse, and diagonal stiffener 
segment are 


Sxsk 

(^xsA:) a ^ / ^ xsk 

Sysk 

(^ysk) a l / ^ysk 

Sxysk 

= (.llyskUKysk 

Si 

= (4,W4 

s 2 

= (4W4 

S 3 

= (4W& 


(16) 


where (e° s *) a /, (4^) a/, ('l xs k)ai and (e£ s *) a f are the allowable membrane strains in the 
skin and stiffeners, respectively. 


2.7 OPTIMIZATION OF GRID-STIFFENED CYLINDERS 


The design variables for a grid-stiffened composite shell are the axial and trans- 
verse stiffener spacings (a, b), the stiffening configuration (ICON), which is the 
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combination of axial, transverse and diagonal stiffeners, the skin laminate ( LAMI ), 
and the height (&), and thickness (t s ) of the stiffener. Except for the height of the 
stiffener, these design variables take on discrete values. 

The genetic algorithm is a method for “evolving” a given design problem to 
a family of near-optimum designs (e.g., see References [10] and [11]). Stochastic 
processes are used to generate an initial population of individual designs and the 
process then applies principles of natural selection and survival of the fittest to find 
improved designs. Furthermore, since the discrete design procedure works with a 
population of designs it can explore a large design space and climb different hills. This 
is a major advantage as the converged solution contains many optima of comparable 
performance. The cost of having a large number of function evaluations is offset by 
the fact that a large number of optimum solutions are now available. The population 
or family of good designs produced by using the genetic algorithm may include the 
global optimal design, rather than a single design. 

2.T.1 Design Problem Definition 

The present design problem is to minimize the weight of a grid-stiffened compos- 
ite circular cylindrical shell given the design loading condition, the length and radius 
of the cylinder, and the material properties for the skin and stiffeners. The design 
variables include stiffener spacings (a, 6), the stacking sequence of the skin, stiffener 
layout, stiffener thickness (f s ), and stiffener height (hi = h 2 = h 3 = h) as shown in 
Figure 1. All stiffeners are assumed to be of the same height and thickness for man- 
ufacturing and assembly reasons. The design sought here is a cylinder of minimum 
weight in a certain design space which buckles globally at the design loads while the 
membrane strains in the skin and the stiffener segments do not exceed the allow- 
able membrane strains (e°, t ) 0 i, (ej,*) 0 i, h° ysk )al and (e° xst ) a i respectively. This design 
problem can be defined by setting up the optimization procedures in the following 
way. First, the global buckling load is assumed to be a scalar multiple of design loads 
and has the form 


(N X , Ny, N X y ) = \g(N\, N2,N\2) 


(17) 
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where Ni, N 2 , and iV 12 , are the applied in-plane prebuckling load. This values rep- 
resent the design loads for the grid-stiffened cylinder. Second, the design constraints 
imposed on panel include 

1. The critical buckling load should be greater than or equal to the design loads, 
that is, A g > 1. 


2. Skin segments should not buckle at the critical buckling load, that is, A > 1. 

3. Stiffener segments should not cripple at the critical buckling load, that is, 

^1? ^3 > T 

4. The membrane strains in the skin segment should be less than or equal (e ° s *.) a /7 

and (q xys k)al that is, S xs k, Syski Sxysk > 1* 

5. The axial membrane strain in the stiffener segment should be greater than or 
equal (e° 5f ) a {, that is, S\, Si, S 3 > 1. 


The general form of each constraint equation is written as 

\ (x- - 1) < o-o 


9, 


(i-l)<0.0 


j = 1 ,...,N C 


(18) 


Finally, the “Fitness” expression based on exterior penalty function approach is 


Fitness = ( 


Q 

F(X,n) 


) = Max 


Q 


W(X) + r,Z»‘ [| ft -(X)|+ ft (X)P 


where X = design variable vector 

F(X, r t ) = Modified objective function 

fK(X) = weight of of cylinder 

n Ef c [|^(X)| + ^(X)] 2 = penalty function 

Q — normalizing constant 

N c — Number of design constraints 

r t - = penalty parameter 

i = generation or iteration cycle in the optimization procedure. 


(19) 
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2.7.2 Design Process Based on Genetic Algorithm 

Implementation of the genetic algorithm is shown schematically in Figure 5. 
The design process begins with a random selection of a specified number of designs 
which comprise the initial population (i.e., first generation) for the genetic algorithm. 
Material properties, radius and length of the cylinder, boundary conditions of the skin 
segment, and design loadings are input to the analysis processor routine. The buckling 
analysis is performed which provides the critical eigenvalues for the global buckling 
response of the grid-stiffened cylinder, the local buckling response of the skin and 
stiffener segments, and the strain level factors of the skin and stiffener segments. The 
weight of the grid-stiffened cylinder is also computed. This procedure is repeated for 
each design configuration in the population. The “fitness” processor then evaluates 
the “fitness” of each design using Equation (19) and assigns a rank based on the fitness 
expression or objective function. The current population of design configurations 
is then processed by the genetic operators (crossover, mutation, and permutation) 
to create a new population of design configurations for the next generations which 
combines the most desirable characteristics of previous generations. Designs from 
previous generations may be replaced by new ones (i.e., children) except for the 
“most fit” designs (i.e., parents) which are always included in the next generation. 
The process is repeated until design convergence is obtained, which is defined herein 
by specifying a maximum number of generations ( NSTOP ) that may occur without 
improvement in the best design. 
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Chapter 3 

USER INSTRUCTIONS 


User instructions for using two FORTRAN codes are provided in this chapter. 
The first code is for analyzing a grid-stiffened cylinder subjected to combined in- 
plane loading. The code provides the global buckling load, local buckling load of skin 
and stiffener segments, and the strain level factors as output. Instructions for using 
this code are given in Section 3.1. The second code is for optimizing a grid-stiffened 
cylinder design subjected to global and local buckling constraints and strength con- 
straints. Instructions for using this code and modifying it to obtain a particular type 
of optimization are given in Section 3.2. 

3.1 ANALYSIS CODE 

The analysis code is found in directory “cylinder/ analysis” and a makefile is used 
to link all the subroutines together. The makefile may have to be modified to account 
for different computer system. Grid-stiffened cylinder with the unit cell geometry as 
shown in Figure 1 can be analyzed by using the code. Skin segments of the grid- 
stiffened cylinder are assumed to be simply-supported. Other boundary conditions of 
the skin segment can be accommodated through simple modifications to the source 
code. The executable for this code is “run” as specified in the makefile. 

3.1.1 Examples for Input and Output file 

An input file for a grid-stiffened cylinder with axial and diagonal stiffeners is 
given. The cylinder is 291.0 in. long, with a radius of 95.5 in., and has an axial 
and transverse stiffener spacings of 8.31428 in. and 14.4588 in., respectively. The 
height and thickness of the stiffener is 0.4125 in. and 0.09 in., respectively. The skin 
laminate has a ply stacking sequence of [±45/02]2s with ply thickness of 0.008 in. The 
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ply orientations are measured from the x-axis in the counter-clockwise directions. The 
stiffener is made of unidirectional material. The material for the skin and stiffener 
is assumed to have the following nominal ply mechanical properties; E\ 1 = 20.2 Msi; 
E 2 2 — 1.9 Msi; G \2 — G\z — G 2 3 = 0.73 Msi and v\ 2 — 0.3. The cylinder is subjected 
to an axial compression loading of N x — 1980 lbs/in. The allowable strain are (e° sjt ) 
— 2428E-06, and = 1092E-06. In this example (e° sA ) — {ll ysk ) = 0.0. The 
input file is named “pan.inp” and is given in List 1. Text after the character “!” are 
comments and need not be included in the actual file. 

Some considerations for the input file are listed below. 

1. The maximum number of plies in a laminate is 50, and the maximum number 
of material is 5. 

2. When a stiffener type (axial, transverse, or diagonal) is not present, its thick- 
ness and ply thickness are entered as zero. But not the material properties 
and its height as shown for the transverse stiffener in the input file (List 1). 

3. When analyzing cylinders stiffened in one axial direction only, one of the stiff- 
ener spacing is redundant. For example, an axially stiffened cylinder will have 
its stiffener spacing specified by the width of the unit cell only. The length of 
the unit cell is entered as the length of the cylinder. Cylinders stiffened in the 
circumferential or transverse directions are not considered herein. 

4. The number of terms or modes (N) used in the analysis is taken from Table 
1 for the local buckling of the skin. The maximum value for N is 100 as 
determined by parameter “nmod” in the file “panel. inc”. Using a value of N 
between 55 and 78 is usually sufficient. The maximum number of terms (N) 
for the global buckling is 25 as determined by parameter “ncyl” in the file 
panel. inc (N = y/ncyl). 

Finally, if the user wishes to change the boundary conditions of the skin segment 
in the analysis, the subroutine “bclocal.f” has to be modified. The output file is 
“pan. out” and is given in List 2. 
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3.2 OPTIMIZATION CODE 

The optimization code will optimize a grid-stiffened cylinder for minimum weight 
subjected to global and local buckling constraints, and strength constraints, and is 
found in directory “cylinder/optimize”. The executable for this code is “run” as 
specified in the makefile. The design variables are 

1. Axial stiffener spacing (a). 

2. Transverse stiffener spacing ( b ). 

3. Stiffener height ( h ). 

4. Stiffener thickness (£ s ). 

5. Skin laminate ( LAMI ). 

6. Stiffening configuration {IGEO). 

The number of design variable is defined by the parameter “n” in the main program 
“main.f’ Each design variable can assume eight discrete values as allowed by the 
FORTRAN code. The eight discrete values of each design variable define the design 
space for optimization. The discrete values for a, 6, h , and t s are supplied through 
the file “inp.gen” which is read by the main program “main.f’. Part of the main 
program where the parameter “n”, the parameter for the population size “m” are 
defined, and the values for a, 6, /i, and t s are read is shown in List 3. The discrete 
values for LAMI and IGEO are given in Table 2. The weight of the cylinder depends 
on the density of the material used. The density of the material, p, is hard-wired in 
subroutine “volume.f’ and can be changed by adjusting the statement “rho = 0.057” 
in the subroutine. 
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Table 2 Design space for design variables ICON and LAM I. 


Integer 

value 

LAMI 

IGEO 

1 

[±45/0] 2s 

axial stiffeners 

2 

[±45/90] 2s 

axial stiffeners’" 

3 

[±45/0/90] 2s 

axial and transverse stiffeners 

4 

(±45/0 2 ] 2iS 

diagonal stiffeners 

5 

[±45/90 2 ] 2s 

axial and diagonal stiffeners 

6 

[±45/0 2 /90] 2s 

transverse and diagonal stiffeners 

7 

[±45/0/90 2 ] 2s 

axial, transverse and diagonal stiffeners 

8 

[±45/0 2 /90 2 ] 2s 

no stiffeners 

* Cylinders with circumferential stiffeners only are not considered. 


The laminate stacking sequence corresponding to various discrete values of 
LAM I are hard-wired in subroutine “cskin.P, and can be changed by modifying 
subroutine “cskin.P. The ply thickness in subroutine “cskin.P is kept constant for all 
laminates, and is read in “cskin.P. Subroutine “cskin.P can be modified by the user 
to accommodate various laminate stacking sequences. The discrete values for I GEO 
are assigned in subroutine “panel. P and part of the code where I GEO are being 
assigned is shown in List 4. Subroutine “geom.P assigns the stiffening configuration 
based on the value of IGEO which is supplied by the main program “main.P. 

Some parameters that may affect the optimization process are 

• The population size “m” is hard-wired in “main.f’. Usually m > 2n, and this 
condition has been found to work well, and m 2n is not recommended ([11]). 

• The probabilities of crossover, mutation, and permutation have been hard- 
wired to 1.0, 0.1, and 0.95 in “main.P ([11]). 

• The termination criteria “NSTOP” is hard- wired in “main.P. The user has 
to experiment with the value of “NSTOP”. Usually the code is run with a 
value of “NSTOP” and then with another value of “NSTOP” greater than 
the previous one. If there is no change in the optimal designs, then the second 
value of “NSTOP” provide a good value as a stopping criteria. For the problem 
under consideration, NSTOP=25, is usually sufficient. 
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• The penalty parameter r* in Equation (19) can either increase at every i th 
generation or can be constant for all generation. In subroutine “panel. f”, r* is 
kept constant at 1000. By commenting the line where “ ainipen = 1000.0”, 
the user can set 

n = 1000 + i 2 (20) 

Keeping r t - constant works very well for the present optimization problem. 

According to Lists 3 and 4, the code has been set up to optimize grid-stiffened 
panel with all design variables active. 

3.2.1 Changing the Type of Optmization 

The user may want to optimize a grid-stiffened cylinder with less number of 
design variables. For example, the skin laminate and the stiffening configuration 
are fixed, and the only design variables are a, 6, /i, and t s . An example of such an 
optimization is provided with the required input files, and the output files from the 
code are explained. 

Consider the cylinder described in Sub-section 3.1.1, the cylinder is to be opti- 
mized for N x = 1980 lbs/in., with design variables being a, 6, h, and t s . The skin 
laminate is [±45/02]2s with a ply thickness of 0.008 in. Only axial and diagonal 
stiffeners are considered and therefore IGEO — 5. The axial stiffener spacing a and 
transverse stiffener spacing b is treated as one design variable i.e., (a, b) is a design 
variable. Values for (a, b) are provided such that the stiffening configuration closely 
approximates an isogrid configuration (i.e., 0 ~ 30°). Hence, there are three design 
variables. In this example, the allowable strains are set to zero, and hence the strength 
constraints are inactive. List 5 and 6 show the appropriate modifications to “main.f” 
and ’’panel.f’ respectively. In List 5, “n” has been changed to 3, and “m” has been 
changed to 8. While in List 6, a “! modify” indicates the line that has been modified. 

The code needs two input files, namely “inp.gen”, and “pan.inp”. The file 
“inp.gen” is read by program “main.f” and it defines the design space for the stiffener 
spacings, and the height and thickness of the stiffener. The file “pan.inp” provides the 


22 


problem parameters for the optimization problem and is read by subroutine “panel.f’ . 
Example for “inp.gen”, and “pan.inp” are given in List 7 and 8 respectively. 

The output files produced by the code are “best. gen” , “on. gen”, and “pan. out”. 
The files “best. gen”, and “on. gen” are produced by program “main.f”, and the file 
“pan. out” is produced by subroutine “panel.f’. The optimal designs ranked according 
to Equation (19) are stored in the file “best.f’ and the convergence history of the 
optimization is stored in the file “on. gen”. The file “best. gen” and “on. gen” for the 
above example is given in Lists 9 and 10 respectively. 

The file “pan. out” stores the information about each design resulting from the 
analysis and is quite large. It is useful in obtaining the buckling loads and other 
information about the optimal designs stored in file “best. gen”. To access information 
about an optimal design, use the value of its fitness (FS) in “best.gen” and locate 
that number (critlb) in “pan. out” using the search option of the unix editor being 
used. For example, the best design, which is the first design in “best.gen” has “FS— 
0.3512084E+00”. Searching for the pattern “0.3512084E-f-00” in “pan. out” will bring 
the cursor to where information about the best design is written. List 11 and 12 give 
information about the first and second optimal designs which have been extracted 
from “pan. out”. 
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LIST OF FILES 


List 1: Example for Input file (pan.inp) 


0.128 ! thickness of skin 

1 ! Number of material 

1 ,20 . 2e6, 1 . 9e6 ,0 .3,0 .73e6 ! Material No., Ell, E22, V12, G12 
16 ! No. of plies 

1.1.0. 008.45.0 ! layer No., Material No., ply thickness, theta 

2.1.0. 008,-45.0 

3.1.0. 008.0.0 

4.1.0. 008.0.0 

5.1.0. 008.45.0 

6.1.0. 008,-45.0 

7.1.0. 008.0.0 

8.1.0. 008.0.0 

9.1.0. 008.0.0 

10.1.0. 008.0.0 

11.1.0. 008,-45.0 

12.1.0. 008.45.0 

13.1.0. 008.0.0 

14.1.0. 008.0.0 

15.1.0. 008,-45.0 

16.1.0. 008.45.0 


♦axial stiffener thickness 
Number of material 
Material No., Ell, E22, V12, G12 
No. of layers 

! layer No., Material No,, ply thickness, theta 
♦transverse stiffener thickness 
Number of material 
Material No., Ell, E22, V12, G12 
No . of layers 

! layer No., Material No., ply thickness, theta 
♦diagonal stiffener thickness 
Number of material 
Material No., Ell, E22, V12, G12 
No. of layers 

1.1.0. 090.0.0 ! layer No., Material No., ply thickness, theta 

0.4125,0.4125,0.4125 ! height of X, Y, D-stiffener 

291.0. 95.5 ! length, radius of cylinder 


0.09 

1 

I,20.2e6,1.9e6,0 . 3, 0 .73e6 
1 

1.1.0. 090.0.0 
0.0 

1 

1 , 20 . 2e6 , 1 . 9e6 ,0 . 3, 0 .73e6 
1 

1 . 1 . 0 . 0 . 0.0 
0.090 

1 

1,20 .2e6,1.9e6, 0.3,0. 73e6 
1 


8.31428,14.4588 

15 

45 

1980.0,0.0,0.0 


length, width, orientation of U. cell 
max m,n in Fourrier series 
# of modes considered for local buckling 
Nx, Ny, Nxy (loading condition) 


2428. 0e-06, 0.0, 0.0, 1092. 0e-06 ! skin_x, y, xy, stiff (allow, strains) 
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List 2: Example for Output file (pan. out) 


SKIN LAMINATE DATA 


STACK THICKNESS 3 0.1280000000000000 
NO. of MATERIAL types 3 1 


m 


[AT. NO. 

El 

E2 

V12 G12 


(psi) 

(psi) 

(psi) 

0 , 

. 2020E+08 

0 . 1900E+07 

0.3000 0 . 7300E+06 

,AYER 

MATERIAL 

THK 

ORIENTATION 

No. 

No. 

(in) 

(deg) 

1 

1 

0.0080 

45.0000 

2 

1 

0.0080 

-45.0000 

3 

1 

0.0080 

0.0000 

4 

1 

0.0080 

0.0000 

5 

1 

0.0080 

45.0000 

6 

1 

0.0080 

-45.0000 

7 

1 

0.0080 

0.0000 

8 

1 

0.0080 

0.0000 

9 

1 

0.0080 

0.0000 

10 

1 

0.0080 

0.0000 

11 

1 

0.0080 

-45.0000 

12 

1 

0.0080 

45.0000 

13 

1 

0.0080 

0.0000 

14 

1 

0.0080 

0.0000 

15 

1 

0.0080 

-45.0000 

16 

1 

0.0080 

45.0000 


X-STIFFENER LAMINATE DATA 


STACK THICKNESS 3 8 . 9999999999999997E-02in 
NO. of MATERIAL types 3 1 


MAT. NO. 


El 
(psi) 

0 . 2020E+08 


E2 

(psi) 


V12 


G12 

(psi) 


0 . 1900E+07 0.3000 0.7300E+06 


LAYER 

No. 

1 


MATERIAL 

No. 

1 


THK 

(in) 

0.0900 


ORIENTATION 

(deg) 

0.0000 


Y-STIFFENER LAMINATE DATA 


STACK THICKNESS 3 0 . OOOOOOOOOOOOOOOOE+OOin 
NO. of MATERIAL types 3 1 


MAT. NO 


El 

(psi) 


E2 

(psi) 


V12 


G12 

(psi) 
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1 0 . 2020E+08 0 . 1900E+07 0.3000 0.7300E+06 

LAYER MATERIAL THK ORIENTATION 

No. No. (in) (deg) 

1 1 0.0000 0.0000 

DIAGONAL-STIFFENER LAMINATE DATA 

STACK THICKNESS- 8 . 9999999999999997E-02in 
NO. of MATERIAL types= 1 

MAT. NO. El E2 V12 G12 

(psi) (psi) (psi) 

1 0 . 2020E+08 0 . 1900E+07 0.3000 0.7300E+06 

LAYER MATERIAL THK ORIENTATION 

No. No. (in) (deg) 

1 1 0.0900 0.0000 

Height of X-stiffener = 0.4125000000000000 

Height of Y-stiffener = 0.4125000000000000 

Height of Dia-stif f ener = 0.4125000000000000 

Length = 291.0000000000000 

Radius = 95.50000000000000 

STIFFENER ORIENTATION (deg) = 29.90030151121361 

UNIT CELL LENGTH (in) = 8.314280000000000 

UNIT CELL WIDTH (in) = 14.45880000000000 

No of MODES FOR CYLINDER = 15 

No of MODES CONSIDERED = 45 

*U* *V* * W * *pX * *pY * 

0000000000 

1010101010 

0101010101 

2020202020 

1111111111 


0808080808 

LOADING MATRIX 

1980.000 0.000 

0.000 0.000 

STRAIN ALLOWABLE 

Max. axial strain in skin = 2 .4280000000000000E-03 

Max. transverse strain in skin = 0 . 0000000000000000E+00 

Max. shear strain in skin — 0 . OOOOOOOOOOOOOOOOE+OO 
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Max. axial strain in stiffener = 1 . 0920000000000001E-03 


END OF INPUT DATA 


SKIN STIFFNESS DATA 


Extensional Stiffness (lbs) 



1725572.065 

365086.081 

0.000 


365086.081 

544372.803 

0.000 


0.000 

0.000 384943.176 


Coupling Stiffness (lbs) 



0.000 

0.000 

0.000 


0.000 

0.000 

0.000 


0.000 

0.000 

0.000 


Bending Stiffness 

(lbs-in) 



1904.344 

647.714 

94.496 


647.714 

896.388 

94.496 


94.496 

94.496 

674.825 


Transverse shear 

stiffness (lbs/in) 



0 . 77867E+05 

0 . 00000E+00 



0 . 00000E+00 

0 . 77867E+05 



STIFFENER EXTENSIONAL COUPLING 

BENDING 

SHEAR 

STIFFNESS 

(lbs) (lbs in) 

(lbs in'‘2) 

(lbs) 


1 0 . 7499E+06 -.2027E+06 0.6540E+05 0.2710E+05 

2 0 . 0000E+00 0 . OOOOE+OO O.OOOOE+OO 0. 0000E+00 

3 0 . 7499E+06 -.2027E+06 0.6540E+05 0.2710E+05 

« ONLY AXIAL & DIAGONAL STIFFENERS » 

Stiffening Parameter (X) = 7 . 0050200635535451E-02 

Stiffening Parameter (Y) - 0 . OOOOOOOOOOOOOOOOE+OO 

Stiffening Parameter (D) = 2 . 89353487754499 19E-02 

cxc = 9 . 2134956352011789E-08 

cxs = 0. OOOOOOOOOOOOOOOOE+OO 

sxc = -3 .4502801522281621E-08 

sxs = 0. OOOOOOOOOOOOOOOOE+OO 

znn = -1 . 8697911092877051E-02 

zstar = -9 . 8630197802043167E-03 

ystar - 0.7139400000000000 

nstep = 10 

cxc = 1.7963803460849354E-07 

cxs = 1.5821236813778613E-07 

sxc = -6.2975256903652235E-07 

sxs = -1.2523730494538431E-07 

znn = -3.6583424817610677E-02 

zstar = -I . 8291712342129864E-02 

ystar » 8 . 2494266793107549E-02 

nstep = 1 


CORRECTED STIFFNESS 


STIFFENER 

EXTENSIONAL 

COUPLING 

BENDING 

SHEAR 

STIFFNESS 

(lbs) 

(lbs in) 

(lbs in~2) 

(lbs) 

1 

0 . 7499E+06 

-. 1968E+06 

0 . 6148E+05 

0 . 2710E+05 

2 

0 . OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

0. OOOOE+OO 

3 

0 . 7499E+06 

1903E+06 

0 . 5824E+05 

0 . 27 10E+05 


EXTENSIONAL SMEARED STIFFNESS 


Stiffeners 

116573.817 

38857.468 

0.000 

stiffeners + skin 
1842145.881 
403943.549 
0.000 


38857.468 

117514.026 

0.000 

403943.549 

661886.829 

0.000 


0.000 

0.000 

38857.468 

0.000 

0.000 

423800.644 


COUPLING SMEARED STIFFNESS 


Stiffeners 

-30481.518 

-9861.939 

0.000 

stiffeners + skin 
-30481.518 
-9861.939 
0.000 


-9861.939 

-29824.800 

0.000 

-9861.939 

-29824.800 

0.000 


0.000 

0.000 

-9861.939 

0.000 

0.000 

-9861.939 


BENDING SMEARED STIFFNESS 


Stiffeners 

9501.363 

3017.777 

0.000 

stiffeners + skin 
11405.708 
3665.491 
94.496 


3017.777 

9126.460 

0.000 

3665.491 

10022.848 

94.496 


0.000 

0.000 

3017.777 

94.496 

94.496 

3692.602 


SMEARED TRANSVERSE SHEAR STIFFNESS (lbs/in) 


stiffeners 

5617.481 

0.000 

stiffener + skin 

5617.481 

0.000 


0.000 

5651.461 

0.000 

5651.461 


BEGIN BUCKLING ANALYSIS 


lcons= 

ierr 


225 

0 
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glam = 0.9713123932539653 

skfx = 1800.514438239891 

skfy = 0.0000000000000000E+00 

skfxy = O.OOOOOOOOOOOOOOOOE+OO 

stdfx = 236.9871816396951 

stdfy = O.OOOOOOOOOOOOOOOOE+OO 

stdfxyx = O.OOOOOOOOOOOOOOOOE+OO 

stdfxyy = O.OOOOOOOOOOOOOOOOE+OO 

stfx = 1913.152055822510 

DETERMINANT = 60.10725583200000 

skilam = 1.248682754551860 

riblx = 1.880641154677955 

ribld = 30.42410048392902 

rlam = 1.000000000000000 

xxl = 1.996763412355778 

yyl = O.OOOOOOOOOOOOOOOOE+OO 

xyl = O.OOOOOOOOOOOOOOOOE+OO 

xstl = 1.037612920905549 

xstl3 = 16.80356955513164 

Volume = 1428.192546820781 ! (weight 
CPU TIME 8.131398 


i Global lambda 
! (Nx)_sk 
! (Ny)_sk 
! (Nxy)_sk 
! N_dx 
! N_dy 
! (N_dxy)_x 
! (N_dxy)_y 
! (N_x) _1 

! lambda.skin 
! lambda_x_stiff 
! lambda.d^stiff 

! S_xsk 
! S_ysk 
! S_xysk 
! S_1 
! S_2 

(lbs) density = 0.057 lbs/in"^ 


o o o o o 
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List 3: Part of main program “main.f” 


C GENETIC ALGORITHM 

C 

C IS(I,J), I being the individual number and J its Jth bits. 
C IS(I } J) bit string number I (old generation) 

C JS(I,J) bit string number I (new generation) 

C CRITLB(I) fitness associated to the individual I 

C FNS(I) normalized fitness of the individual I 

C M population size 

C NLA maximum number of layers 

C N number of bits in a string (=NLA/4) 

C PC probability of implementing crossover 

C PM probability of implementing mutation 

C PP probability of permutation 

C PRI probability of inversion 

C ITER iteration (generation) number 
C 
C 

parameter (n=6 ,m=12 ,nn=m*2) 
c 

c n = number of design variables 

c m = number of design in each group 

c nn= just for dimension ( parameter) 
c 

DIMENSION IS(nn,n), FNS(IOO), RD(n) ,JS(m,n), 

&CRITLB (100) ,CRITZ(100) ,aas(8) ,bbs(8) ,hhl(8) ,tthkl(8) , 
&ISK(50,m),FSK(50) ,NGEN(50) ,alength(n) 
real*4 tcp2(2) 
c 

COMMON/PIE/PI 
C0MM0N/MATGE0/T , NLA 
common /function/critlb 
c 

OPEN (UNIT=15, FILE- 1 on. gen’ ) 

OPEN (UNIT-12 , FILE='best.gen>) 

OPEN (UNIT-9 , FILE='inp.gen') 

C 

c nla=16 

c N=NLA/4 

n-15 


Maximum number of generations 

read (9,*) (aas(ij) , i j = 1,8) 
read (9,*) (bbs(ij) , ij-1 ,8) 
read (9,*) (hhl(ij) ,ij=l,8) 
read (9,*) (tthkl(ij) , i j = 1 ,8) 

C 

LTT=300 

C 


o o o o o o o ooooo oooon o o 
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Genetic parameters 

PC-1 . 00D+00 
PM=0 . 10D+00 
PP=0. 95+00 
M-16 


NFT is the number of evaluations of the objective function 
without improvement before the search stops. 

NFT=150 


Initialization of the stopping criterion 

NCRIT=0 
0PTI=0 . D+00 

NSTQP is the maximum number of generations without 
improvement . 

NST0P=20 

Initialization of the parameters of the subroutine STORE 
before the first call. 


call dtime(tcp2) 

write(12 ,*) 'CPU TIME =\tcp2(l) 

CL0SE(12) 

END 
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List 4: Part of subroutine “panel.f” 

subroutine panel (io , is ,critlb , ainipen,nn,nd,aas } bbs ,hhl , 
& tthkl) 


include "opt.inc" 
include "panel. inc" 
integer is(nn,nd) 
c 

c --------------- 

c UNIT 5 FOR THE INPUT DATA FILE 


open (5 ,f ile= ’pan . inp ’ ) 
rewind 5 


UNIT 6 IS FOR THE OUTPUT DATA FILE 
open (6 ,f ile=’pan.out ’ ) 


c 

c 

c 

c 

c 

c 

c 


c 


clen = aas(is(io , 1) ) 
cwid ~ bbs (is(io , 2) ) 
hi = hhl(is(io,3)) 
h2 = hi 
h3 = hi 

thick = tthkl(is(io,4)) 
igeo =* is(io, 
laini = is(io,5) 
igeo = is(io,6) 

call geom(sthl , sth2 , sth3 , igeo , thick) 

write (6 ,*) > 

write (6 ,*) ; Laminate salami 

write (6 ,*) ’Stiffener thickness = ’ ,thick 

write (6,*) ’Stiffener height = ’ ,hl 


[1] READING ALL INPUT DATA FOR A LAMINATE 

USING SUBROUTINE ISKIN 


call iskin(sthk,nmsk,exsk,eysk,vsk,gsk,amatsk,nsk, 
& lnosk,msk,tsk,thesk) 

call cskin ( sthk , nmsk , exsk , eysk , vsk , gsk , amat sk , nsk , 
& lnosk,msk,tsk,thesk,lami) 


write (6 ,51)critlb(io) 
f ormat ( ’ critlb - ’el4.7) 

write (6 , *) ’ 

return 

end 


51 
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List 5: Modifications to main program “main.f” 


C GENETIC ALGORITHM 

C 

C IS(I,J), I being the individual number and J its Jth bits. 

C IS(I,J) bit string number I (old generation) 

C JS(I,J) bit string number I (new generation) 

C CRITLB(I) fitness associated to the individual I 
C FNS(I) normalized fitness of the individual I 
C M population size 

C NLA maximum number of layers 

C N number of bits in a string (-NLA/4) 

C PC probability of implementing crossover 

C PM probability of implementing mutation 

C PP probability of permutation 

C PRI probability of inversion 

C ITER iteration (generation) number 
C 

O* ****************************************************** *********** 
parameter (n=3 ,m=8 ,nn=m*2) 
c 

c n = number of design variables 

c m = number of design in each group 

c nn= just for dimension ( parameter) 
c 

DIMENSION IS (nn,n) , FNS(lOO), RD(n) ,JS(m,n), 

&CRITLB(100) ,CRITZ(100) ,aas(8) ,bbs(8) ,hhl(8) ,tthkl(8) , 

&ISK(50 ,m) ,FSK(50) ,NGEN(50) ,alength(n) 
real*4 tcp2(2) 
c 

COMMON/PIE/PI 
COMMON/MATGEO/T , NLA 
common /function/critlb 
c 

OPEN (UNIT-15 , FILE= , on.genO 
OPEN (UNIT=12 , FILE= ' best .gen ; ) 

OPEN (UNIT=9 , FILE= J inp.gen') 


read (9,*) (aas(ij) ,ij=l ,8) 
read (9,*) (bbs(ij) ,ij=l ,8) 
read (9,*) (hhl(ij ) ,ij=l ,8) 
read (9,*) (tthkl(ij) ,ij=l ,8) 


END 
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List 6: Modifications to subroutine “panel. f” 

subroutine panel (io , is , critlb, ainipen,nn,nd, aas ,bbs ,hhl , 
& tthkl) 


include "opt.inc" 
include "panel. inc" 
integer is(nn,nd) 
c 

c UNIT 5 FOR THE INPUT DATA FILE 

c 

open (5 , f ile= ' pan . inp ' ) 
rewind 5 
c 
c 

c UNIT 6 IS FOR THE OUTPUT DATA FILE 

c 

open (6 , f ile= ' pan . out ' ) 


c 


c 

c 

c 

c 

c 

c 

c 


c 


clen = aas(is(io, 1) ) ! modify 
cwid = bbs(is(io, 1) ) ! modify 
igeo = 5 ! modify 
lami - 4 ! modify 
hi = hhl(is(io,2)) ! modify 
h2 - hi ! modify 
h3 = hi ! modify 
thick = tthkl(is(io,3)) ! modify 


call geom(sthl , sth2 , sth3 , igeo , thick) 

write (6,*) J 

write(6,*) 'Laminate =, ,lami 
write(6,*) 'Stiffener thickness = ', thick 
write(6,*) 'Stiff ener height = ' ,hl 


[1] READING ALL INPUT DATA FOR A LAMINATE 

USING SUBROUTINE ISKIN 


call iskin ( st hk , nmsk , exsk , ey sk , vsk , gsk , amat sk , nsk , 
& lnosk,msk,tsk,thesk) 

call cskin ( st hk, nmsk, exsk, eysk, vsk, gsk, amat sk, nsk. 
Sc lnosk,msk,tsk,thesk,lami) 


END 



List 7: Input file for optimization (inp.gen). 


7 . 4615385 , 7 . 6578947 , 7 . 86486 , 8 . 08333 , 8 .31428 ,8.559, 

8.818181,9.09375 ! a 

12 . 904760 , 13 . 334316 , 13 . 637386 , 13 . 954510 , 14.458800 , 

14.815906,15.385749,15.790637 ! b 

0.40,0.4125,0.425,0.4375,0.45,0.4625,0.475,0.4875 ! h 

0.048,0.054,0.06,0.066,0.072,0.078,0.084,0.090 ! t 


List 8: Input file for problem parameters (pan.inp). 


1 

1,20. 2e6 , 1 . 9e6 ,0 
0.008 
1 

1 ,20.2e6, 1 . 9e6 ,0 
1 

1 , 1 , 0.0 
1 

1,20. 2e6 , 1 . 9e6 ,0 
1 

1 , 1 , 0.0 
1 

I,20.2e6,1.9e6,0 

1 

1 , 1 , 0.0 

291.0. 95.5 
15 

45 

1980.0. 0.0.0.0.0 
0 . 0 , 0 . 0 , 0 . 0 , 0.0 


Number of material 

3 , 0 .73e6 ! Material No., Ell, E22, V12, G12 
ply thickness for skin laminate 
♦Number of material (X-st if f ener) 

3,0.73e6 ! Material No., Ell, E22 , V12, G12 
Number of plies 

layer No., material No., theta 
♦Number of material (Y-stif f ener) 

3,0 .73e6 ! Material No., Ell, E22 , V12, G12 
Number of plies 

layer No., material No., theta 
♦Number of material (D-st iff ener) 

3,0 .73e6 ! Material No., Ell, E22, V12, G12 
Number of plies 

layer No., material No., theta 
length, radius of cylinder 
max m,n in Fourrier series 
# of modes considered for local buckling 
Nx, Ny, Nxy 

skin_x, _y, _xy, stiff_x (strain allowable) 
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List 9: Output file containing optimal designs (best. gen). 


POPULATION SIZE 3 8 CROSSOVER PROB. =1.000 
MUTATION PROB . =0 . 100 PERMUTATION PROB . =0 . 950 
BESTS DESIGNS AFTER 246 EVALUATIONS OF THE OF 
FS= 0 . 3512084E+00 GENERATION 3 15 

7 3 8 

FS= 0 . 3500683E+00 GENERATION 3 8 

5 2 8 

FS= 0 . 3499497E+00 GENERATION 3 13 

6 3 8 

FS= 0 . 3494749E+00 GENERATION 3 9 

8 3 8 

FS= 0.3489305E+00 GENERATION 3 4 
5 3 8 

FS= 0 . 3478998E+00 GENERATION 3 34 
8 5 6 

FS= 0 . 3478000E+00 GENERATION 3 11 

5 4 8 

FS= 0 . 3477422E+00 GENERATION 3 33 

6 5 8 

FS= 0.3456823E+00 GENERATION 3 1 
2 3 8 

FS= 0 . 3455608E+00 GENERATION 3 7 


CPU TIME 3 1985.198 
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List 10: Onput file containing convergence history (best. gen). 


population size= 8 crossover prob. =1.000 
mutation prob. =0.100 permutation prob. =0.950 
Iteration Average Best 

1 0 . 1474611E+00 0 . 3456823E+00 

2 0 . 1024948E+00 0 . 3456823E+00 

3 0 . 1993546E+00 0 . 3456823E+00 

4 0 . 2728501E+00 0 . 3489305E+00 

5 0 . 2851492E+00 0 . 3489305E+00 

6 0 . 3256674E+00 0 . 3489305E+00 

7 0 . 2783510E+00 0 . 3489305E+00 

8 0 . 2561589E+00 0 . 3500683E+00 

9 0 . 30331 18E+00 0 . 3500683E+00 

10 0 . 2970200E+00 0 . 3500683E+00 

11 0 . 3478100E+00 0 . 3500683E+00 

12 0 . 3478156E+00 0 . 3500683E+00 

13 0 . 262767 1E+00 0 . 3500683E+00 

14 0 . 2576738E+00 0 . 3500683E+00 

15 0 . 3180880E+00 0 . 3512084E+00 

16 0 . 3490608E+00 0 . 3512084E+00 

17 0 . 3021455E+00 0 . 3512084E+00 

18 0 . 3403848E+00 0 . 3512084E+00 

19 0 . 3053899E+00 0 . 3512084E+00 

20 0 . 2841141E+00 0 . 3512084E+00 

21 0 . 3408947E+00 0 . 3512084E+00 

22 0 . 3505790E+00 0 . 3512084E+00 

23 0 . 3286393E+00 0 . 3512084E+00 

24 0 . 3287967E+00 0 . 3512084E+00 

25 0 . 3507364E+00 0 . 3512084E+00 

26 0.3414941E+00 0 . 3512084E+00 

27 0 . 3363042E+00 0 . 3512084E+00 

28 0 . 33127 16E+00 0 . 3512084E+00 

29 0 . 3270619E+00 0 . 3512084E+00 

30 0.3008105E+00 0 . 3512084E+00 

31 0.3288947E+00 0 . 3512084E+00 

32 0 . 3504603E+00 0 . 3512084E+00 

33 0 . 3406868E+00 0 . 3512084E+00 

34 0 . 3283831E+00 0 .3512084E+00 

FINAL POPULATION AFTER 246 EVALUATIONS OF THE O.F. 
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List 11: Information about first optimal design from (pan.out). 


Laminate = 4 

Stiffener thickness = 8 . 9999999999999997E-02 

Stiffener height = 0.4250000000000000 

STIFFENER ORIENTATION (deg) = 29.81872703828845 

UNIT CELL LENGTH (in) = 8.818180999999999 

UNIT CELL WIDTH (in) * 15.38574900000000 

« ONLY AXIAL & DIAGONAL STIFFENERS » 

Stiffening Parameter (X) = 6 . 7824713502575448E-02 

Stiffening Parameter (Y) = 0 . OOOOOOOOOOOOOOOOE+OO 

Stiffening Parameter (D) = 2 . 7831146619772991E-02 

IFLAG = 5 

znn = -1 . 8748481539636742E-02 

zstar = -9.8976731692590435E-03 

ystar = 0.7602874500000001 

nstep - 10 

znn = -3.6390656669926545E-02 

zstar « -1.8195328248661609E-02 

ystar = 8 . 7768143721372954E-02 

nstep = 1 

lcons= 225 

ierr 0 

Global lambda = 1.004465200334541 

penalty factor = 100000.0000000000 

skfx = 1865.845124687763 

skfy = 0. 0000000000000000E+00 

skfxy = 0. OOOOOOOOOOOOOOOOE+OO 

stdfx = 243.7657369767708 

stdfy = 0. OOOOOOOOOOOOOOOOE+OO 

stdfxyx = 0. OOOOOOOOOOOOOOOOE+OO 

stdfxyy * 0 . OOOOOOOOOOOOOOOOE+OO 

stf x = 1982.569736920487 

DETERMINANT = 67.83715975128450 

skilam = 1.075033581289790 

riblx = 1.698565494963615 

ribld = 27.73106058929319 

rlam = 1.000000000000000 

xxl a 0. OOOOOOOOOOOOOOOOE+OO 

yyl = 0 . OOOOOOOOOOOOOOOOE+OO 

xyl = 0. OOOOOOOOOOOOOOOOE+OO 

xstll a 0. OOOOOOOOOOOOOOOOE+OO 

xstl3 = 0. OOOOOOOOOOOOOOOOE+OO 

Volume = 1423.555983291840 

critlb = 0.3512084E+00 
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List 12: Information about second optimal design from (pan. out). 


Laminate = 4 

Stiffener thickness = 8 . 9999999999999997E-02 

Stiffener height - 0.4125000000000000 

STIFFENER ORIENTATION (deg) = 29.90030151121361 

UNIT CELL LENGTH (in) = 8.314280000000000 

UNIT CELL WIDTH (in) = 14.45880000000000 

« ONLY AXIAL & DIAGONAL STIFFENERS » 

Stiffening Parameter (X) = 7 . 005020063553545 IE-02 

Stiffening Parameter (Y) = 0 . OOOOOOOOOOOOOOOOE+OO 

Stiffening Parameter (D) = 2 . 8935348775449919E-02 

IFLAG = 5 

znn = -1 . 8697911092877051E-02 

zstar = -9 . 8630197802043167E-03 

ystar = 0.7139400000000000 

nstep = 10 

znn = -3 . 6583424817610677E-02 

zstar = -1.8291712342129864E-02 

ystar = 8 . 2494266793107549E-02 

nstep = 1 

lcons= 225 

ierr 0 

Global lambda * 1.001613221328421 

penalty factor = 100000.0000000000 

skfx = 1856.682854104444 

skfy - 0. 0000000000000000E+00 

skfxy = 0. OOOOOOOOOOOOOOOOE+OO 

stdfx = 244.3801768249596 

stdfy = 0. OOOOOOOOOOOOOOOOE+OO 

stdfxyx = 0. OOOOOOOOOOOOOOOOE+OO 

stdfxyy * 0 . OOOOOOOOOOOOOOOOE+OO 

stf x = 1972.834287745410 

DETERMINANT = 60.10725583200000 

skilam = 1.210907572842309 

riblx = 1.823747951708783 

ribld = 29.50370983966332 

rlam = 1.000000000000000 

xxl = 0. OOOOOOOOOOOOOOOOE+OO 

yyl = 0. OOOOOOOOOOOOOOOOE+OO 

xyl = 0. OOOOOOOOOOOOOOOOE+OO 

xstll = 0. OOOOOOOOOOOOOOOOE+OO 

xstl3 = 0. OOOOOOOOOOOOOOOOE+OO 

Volume = 1428.192546820781 

critlb = 0 . 3500683E+00 
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Figure 2 Semi-infinite plate model for skin-stiffener element. 




Figure 3 Typical profile for skin-stiffener element neutral surface. 
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Figure 4 Physical and computational domain for plate geometries. 




Figure 5 Flow chart showing the design optimization process. 




APPENDIX C 


Documentation for the Analysis Code 
for the Buckling of Variable-Curvature Anisotropic Panels 


• Theoretical Foundations 

• User Instructions 

• Sample Problems 


15 


BUCKLING ANALYSIS OF 
ANISOTROPIC 
VARIABLE-CURVATURE 
PANELS AND SHELLS 


Theoretical Foundations and 
User’s instructions for a FORTRAN Code 



Navin Jaunky 

Department of Aerospace Engineering 
Old Dominion University 
Norfolk, VA 23529-0247 


1 


TABLE OF CONTENTS 

TABLE OF CONTENTS i 

LIST OF FIGURES ii 

LIST OF FILES iii 

LIST OF TABLES iv 

NOMENCLATURE v 

Chapter page 

1 INTRODUCTION 1 

2 THEORETICAL FOUNDATIONS 2 

2.1 GEOMETRY OF VARIABLE CURVATURE PANEL 3 

2.2 BEZIER POLYNOMIALS 3 

2.3 CONTINUITIES ALONG SEGMENT JUNCTURES 4 

2.4 MINIMUM POTENTIAL ENERGY PRINCIPLE 6 

3 USER INSTRUCTIONS 8 

3.1 APPLICATION OF BOUNDARY CONDITIONS 8 

3.2 CYLINDRICAL PANEL 9 

3.3 WING LEADING-EDGE PANEL ‘ 10 

3.4 NON-CIRCULAR FUSELAGE 11 

REFERENCES 19 


i 



LIST OF FIGURES 


1 Coordinate system and geometry of shell with variable curvature 22 

2 Control points for joining shell segments 22 

3 Geometry and boundary conditions of curved composite panel 23 

4 Geometry, dimensions and boundary conditions of composite wing 

leading-edge panel 23 

5 Boundary conditions for quarter symmetric model of composite wing 

leading-edge panel 24 

6 Buckling load interaction curve for the composite wing-leading edge panel. . . 24 

7 Geometry of non-circular fuselage 25 

8 Boundary conditions for quarter symmetric model of non-circular fuselage. . . 25 


n 


LIST OF FILES 


1 Input files for cylindrical panel 13 

2 Input file for wing leading-edge panel (full model) 15 

3 Input file for wing leading-edge panel (quarter symmetric model) 17 

4 Input file for non-circular fuselage (isotropic) 18 



LIST OF TABLES 


1 Size of stiffness matrices for panel and shell for increasing 

number of segments 20 

2 Comparison of buckling loads results for composite curved panel 20 

3 Comparisons of buckling loads results for wing leading-edge panel 21 

4 Comparisons of buckling loads results for non-circular fuselage 21 


IV 



Nomenclature 


u 0, Vo 
w 

4*xi fiy 
x, y , 2 

fi(n,v) 

Qij 

t, V 
L 

Si 


displacement along axial and transverse directions 

displacement along radial direction 

rotations of normals to middle surface or curvatures 

axial, circumferential, and normal coordinates 

Bezier polynomial of order n 

Bezier control points 

non-dimensional coordinates 

length of segment 

width of segment I 


R 


radius of curvature 
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ANALYSIS AND DESIGN OF FUSELAGE STRUCTURES 
INCLUDING RESIDUAL STRENGTH PREDICTION METHODOLOGY 

Final Report 

for NASA Grant NAG-1-1588 

by 

Norman F. Knight, Jr. 

Department of Aerospace Engineering 
College of Engineering and Technology 
Old Dominion University 

SUMMARY 

The goal of this research project is to develop and assess methodologies for the design 
and analysis of fuselage structures accounting for residual strength. Two primary objectives 
are included in this research activity: development of structural analysis methodology 
for predicting residual strength of fuselage shell-type structures; and the development 
of accurate, efficient analysis, design and optimization tool for fuselage shell structures. 
Assessment of these tools for robustness, efficient, and usage in a fuselage shell design 
environment will be integrated with these two primary research objectives. 

This research activity extended over a period of four years from January 1, 1994 un- 
til December 31, 1997 with a no-cost extension granted until March 31, 1998. Over the 
course of the grant, the principal investigator, graduate students, a post doctoral research 
assistant, and two research scientists were supported at various levels and during different 
time periods. This research produced nine conference papers, seven archival journal pa- 
pers, one NASA report, one Ph. D. dissertation and four Master’s theses. The research 
computer codes developed under this grant have been documented (see appendices), and 
a distribution CD with FORTRAN source code, sample problems, and postscript versions 
of the documentation is provided. 

BACKGROUND 

The design of aerospace structures generally results in light-weight structural designs, 
advanced structural materials and fabrication concepts, and highly stressed systems. The 
goals of aerospace structural design are to meet the design requirements for the operating 
conditions and flight envelope of the vehicle, adequate service life and damage tolerance, 
and reasonable manufacturing cost. Aerospace fuselage structures are generally subject to 
internal pressure loadings, thermal cycling, bending, axial, and shear loadings, and fatigue 
over its intended life cycle. Fuselage structures involve flat and curved stiffened and un- 
stiffened panels with and without cutouts that are interconnected by frames, stringers and 
bulkheads. Damage tolerance issues associated with fuselage structures have been stud- 
ied by several researchers. The two most common in fuselage structures are longitudinal 
cracks under hoop stresses induced by internal pressure loading and circumferential cracks 
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intersections are only coarsely approximately. Discrete stiffener modeling of these panels 
provided approximately the same level of detail in the response prediction at a fraction of 
the modeling and computational costs. Use of the smeared stiffener theory in STAGS was 
also consider; however, this formulation does not account for any skin-stiffener interaction 
- it simply modifies the skin stiffness by uniformly “smearing” the stiffener stiffness across 
the skin and accounts for eccentricity of the stiffeners. An improved theory is needed for 
preliminary design and analysis tools. 

Shell Theories. As part of our studies on the design of cylindrical panels, various shell 
theories were examined including the Sanders-Koiter theory, the Love theory, and the Don- 
nell theory (see Dl). All are implemented using “tracer” coefficients in the analysis. In 
nearly all cases considered, the buckling results from each theory were in very good agree- 
ment. However, for angle-ply laminates with increasing anisotropy and certain winding 
angles, the results predicted by the Donnell theory were significantly different than those 
obtained using the Sanders-Koiter or Love theory. This was confirmed by STAGS finite 
element analyses. As the winding angle changes, changes in the buckling mode shape occur 
for which Donnell’s theory is not accurate. As we move towards the automated design op- 
timization process, these changes in behavior and their impact on the spatial discretization 
requirements needs to be understood. Results from these investigations are reported in a 
journal paper that has recently been accepted for publication (see JP6). 

Further studies on the influence of which shell theory to use for buckling of cylindrical 
shells were conducted using the PANDA2 computer code from Dr. David Bushnell of the 
Lockheed-Martin Advanced Technology Center. These studies verified the previous results 
for axially compressed cylinders. In addition other loading conditions were considered: 
external pressure, in-plane shear loading and hydrostatic loading. In each of these loading 
cases, Donnell’s theory and Sander’s theory were found to be in good agreement for all 
values of the fiber winding angle considered (see T3). 

Variational Formulation. The original plan was to develop a variational formulation 
of a damaged structure; however, it became necessary to develop a better understanding 
of other local discontinuities and their influence on local stress distributions which then 
contribute to the prebuckling stress state. Close examination of the skin-stiffener inter- 
section region revealed that the traditional assumptions used in smeared stiffener theory 
for preliminary design and sizing of stiffened panels constrained for global overall buckling 
did not account for any interaction between the skin and the stiffener. This local stiffness 
discontinuity results in a shift of the neutral surface of the stiffened panel away from the 
skin middle surface as the stiffener is approached. As a result, the load redistribution 
between the skin and the stiffener is changed. The formulation and its comparison with 
traditional smeared stiffener theory is given in four publications (see CP3, JP3, R1 and 
Dl). 

Progressive Failure Analysis. This activity formed the basis of another Master’s thesis 
(Mr. David W. Sleight, a NASA employee) under the direction of Dr. Knight (see Tl). 
The progressive failure analysis methodology involved the use of COMET and its generic 
constitutive processor (GCP) for the evaluation of point stress failure criteria, degradation 
of ply- level material properties, re-calculation of new laminate stiffnesses, and saving the 
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path-dependent historical material data. The methodology employed various solution pro- 
cedures and processors from COMET which have functionally equivalent routines in the 
current version of STAGS. 

The basic steps of the progressive failure methodology includes the following steps. 
First, a nonlinear solution is obtained at a given load step while holding the material data 
constant during the nonlinear iteration. Then the stress recovery step is performed using 
the element stress resultants (these should be the best available since they are provided 
by the element developer). Next using these stress resultants, the midplane strains and 
changes in curvature are computed and used to calculate the point strains and stresses 
through the thickness of the laminate - only in-plane components are used. Given these 
point values, various failure criteria can be assessed including maximum strain, Chris- 
tensen’s criteria, and Hashin’s criteria. Next, the ply discounting method is used to de- 
grade lamina material data and new laminate stiffness coefficients are computed. Historical 
data are saved and the next iteration begins. 

Several standard test cases with available experimental data were used to verify the 
methodology and its implementation. These included the rail-shear problem and the 
tension-loaded open hole problem. Then compression loaded panels were analyzed and 
the use of these failure models on failure load prediction was performed. Overall good per- 
formance of the progressive failure analysis methodology was obtained as reported in the 
1997 conference paper (see CP5). Sensitivity of the results to material allowable values 
is also noted. In most cases nominal values were used since experimentally determined 
values were unavailable. 

Migration of this methodology to the STAGS finite element code should be possible 
provided that the STAGS code has the GCP features available in COMET. There are 
three critical aspects of this type of analysis. Two are mechanics related: failure mode 
representation and detection and material degradation modeling. The third aspect is the 
organization and preservation of the path-dependent historical data (similar to an elasto- 
plastic analysis). The details of the COMET implementation are in Mr. Sleight’s Master’s 
thesis. 

Shell Analysis for Design and Optimization 

Shell analysis for design and optimization at the preliminary design stage must con- 
sider trade-offs between computational effort and accuracy. The analysis tasks are embed- 
ded within the design optimization iteration process and are frequently performed tens or 
hundreds or even thousands of times for a given design problem. Even though the com- 
puter systems are significantly more powerful today than even five years ago, there still 
are insufficient computational resources to embed detailed finite element structural models 
within the design optimization loop. This limitation is due in part to the fact that often 
times the design optimization procedure requires a geometry change which then requires 
a finite element model change which then requires an engineer in the loop. 

To this end, robust and cost-effective analysis methods are being developed and inte- 
grated into design optimization programs. This research grant considers three independent 
strategies. The first strategy is the use of VICON (or VICON-OPT) which is based on 
extensions to the early work in PASCO and VIPASA. It features exact stiffness matrices 
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for the structural model and has a very accurate eigenvalue solver. This strategy is appli- 
cable to prismatic structures in general. The second strategy is to develop a new family 
of analysis and design tools for general grid-stiffened composite plate and shell structures. 
Significant amount of effort is devoted to this strategy. The third strategy is to explore 
the use of PANDA2 for the design and optimization of sandwich plates and shells account- 
ing for local failures in the sandwich structure as well as possibly postbuckling strength. 
Various aspects of these three strategies are summarized next. 

VICON Program. The VICON computer program for stress and buckling analysis of 
composite panels has been under development for a number of years. It makes use of exact 
stiffness matrices that produce accurate results for any prismatic assembly of flat plates 
without any user requirement to generate a finite element grid. It is also computationally 
efficient achieving results as much as an order of magnitude faster than conventional general 
purpose finite element codes. This efficiency has allowed the development of a design 
capability where the dimensions of a panel (plate breadths, layer thicknesses and ply 
angles) may be adjusted to achieve a minimum mass panel. 

The VICON computer program for analysis and design of plate assemblies has been 
improved in both capability and efficiency. A six-week period was spent by Dr. Ander- 
son at the University of Wales working with the co-developers of the program. The new 
capabilities that have been developed under this grant were combined with the new capa- 
bilities developed at the University of Wales to result in one program having the combined 
capabilities developed at both sites. In discussions with Professor Fred Williams, devel- 
oped a new approach to sensitivity calculations that should result in increased accuracy 
and faster solution times than current method. This was implemented in the program 
and evaluated while at the University of Wales. Other efficiency improvements developed 
under the grant were given extensive evaluation and checks and a number of bugs fixed 
to insure compatibility with new features of the program developed at the University of 
Wales. 

Treatment of Curved Plates. A number of different shell theories have been examined 
to determine the feasibility of incorporating an exact stiffness matrix for a curved plate 
into the VICON program. Use will be made of a numerical method developed for flat 
plates having transverse shear deformation that is in the present program. It is necessary 
that the equations have a certain format for the existing method to be directly applicable. 
It has been found that the theories based on either physical strains or tensor strains and 
neglecting inplane shear and transverse loadings for the inplane equilibrium will satisfy 
this format requirement. The present program which treats only flat plates exactly is 
based on tensor strains with the same neglect of the inplane shear and transverse loadings 
on inplane equilibrium. The theories for the curved plate can be implemented by simple 
modifications of a 10 x 10 matrix in the existing program. The work was accomplished by 
a Master’s graduate student (Mr. David McGowan, a NASA employee) under the direction 
of Dr. M. S. Anderson. He defended his thesis during the Spring 1997 semester (see T2) 
and co-authored a 1997 conference paper on it as well (see CP6). 

Effect of Axial Load Application Point, The equations necessary for the calculation 
of the additional bending that occurs for this case have been developed and implemented 
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in the program. The location of the load application point can be prescribed in input data 
and through linking equations can be made a function of the dimensions of the panel so 
that changes that might occur during the design process can be included. 

Panels having Postbuckled Strength. The original plan was to try to use the program 
P BUCKLE to account for post buckling strength of a panel structure. In looking at the 
capabilities of PBUCKLE it appeared to be limited to certain specific cross-sections rather 
than the general capability of VICON. An alternate idea has been implemented that retains 
the generality of VICON and has the following features: 

1. The designer can choose a load level, less than ultimate, at which buckling will 
occur. 

2. The designer can choose which part of the panel will undergo significant buckling. 

3. The method is quite simple and executes in the same time as would be required 
to design a buckle free panel. 

4. The method requires an assumption of the reduction of the inplane stiffnesses of 
the plates that have been selected to buckle at less than ultimate load. A good 
estimate of this reduction can be obtained from published results on the post 
buckling response of individual plates. 

The method is just reaching its final stages of development and has been applied to 
the design of a metal zee-stiffened panel having a buckling load two thirds of its ultimate 
load. The optimized mass of this panel was found to be about 12% less than if it were 
required to be a non-buckling design at ultimate load. The panel was analyzed with the 
STAGS program to check the validity of the simplifying assumptions made in the VICON 
analysis. The ultimate load determined by the STAGS analysis was 10% higher than the 
design load which is encouraging that the proposed method might be a valid way to design 
panels having postbuckling strength. These results were reported in a 1997 conference 
paper by Dr. Anderson (see CP7). 

Grid-Stiffened Composite Panels. A collection of analysis tools has been developed 
which include a Rayleigh- Ritz approach for linear buckling analysis of the overall panel, 
Rayleigh- Ritz approach for linear buckling of the skin segments locally, and stiffener crip- 
pling assessment. As part of the design process, a global buckling analysis capability was 
developed. Based on these global buckling load levels, the skin segments and stiffener 
segments were individually analyzed to determine if their “local” buckling load exceeded 
the buckling load determined at the global level. If so, then this configuration was an ac- 
ceptable design and its weight was computed. If not, then this configuration was penalized 
as an unacceptable design. 

These analysis methods account for anisotropic material behavior, transverse shear 
deformation effects, and skin-stiffener interaction at the global level through an improved 
smeared stiffener theory. Local buckling analyses were developed to handle different plan- 
form shapes of the skin segments between stiffeners. These methods address general quadri- 
lateral and triangular shaped skin segments and have been verified with classical solutions 
and finite element calculations. Two conference papers (CPI, CP2) and two journal papers 
(JP1, JP2) reported this work. Special attention was also given to the smeared stiffener 
modeling theory used in the global buckling analysis. A new improved theory was devel- 
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oped and implemented. A 1995 conference paper (CP3), 1996 journal paper (JP3), and 
1995 NASA report (Rl) resulted from this effort. 

These analysis tools are integrated together and combined with a genetic algorithm 
for the “evolution” of a family of best designs. The implementation of the genetic algo- 
rithm used in this research is based on the software obtained from Prof. R. Haftka of the 
University of Florida and is gratefully acknowledged. The genetic algorithm or GA is ideal 
for this class of structures because of the inherent discrete nature of the design variables. 
That is, discrete choices exist for stiffener pattern (axial, transverse, diagonal, or selected 
combinations), skin lamination pattern, material type, and stiffener spacings. Integration 
of these various analysis methods with a design strategy based on the genetic algorithm 
was developed and documented (see Appendix A). These results were reported in a 1996 
conference paper (CP4) and a 1998 journal paper (JP4). This research effort formed the 
basis of Mr. Navin Jaunky’s Ph. D. dissertation under Dr. Knight’s direction and was 
completed in December 1995 (see Dl). Dr. Jaunky has continued with this research grant 
since that time as a post-doctoral research associate. 

The source code, sample problems, and documentation (same as Appendix A) are 
available on a compact disk (CD). 

Grid-Stiffened Composite Circular Cylinders. This research effort is focussed on the 
optimal design of general stiffened circular cylinders. In addition to the global buckling 
constraints, an exploratory study has been performed to determine the effect of strength 
constraints in finding an optimal design. Strain allowables are incorporated in to the anal- 
ysis and compared with the computed strains in the stiffeners and panel skins for each 
design configuration during each generation of the genetic algorithm execution. Prelim- 
inary results indicate that the optimal design configuration without strength constraints 
generally leads to a grid- stiffened cylinder which in itself is very redundant in its load 
paths. Design configurations with strength constraints are essentially the same as those 
without the strength constraints except for a slight weight penalty for the case of only 
axial loads. For the combined load cases (e.g., hoop loads or torsion), this is not the case. 
Integration of these various analysis methods with a design strategy based on the genetic 
algorithm was developed and documented (see Appendix B). These results were reported 
in a 1997 conference paper (CP8), and a journal paper (galley proofs have been reviewed, 
publication is pending, JP5). 

The source code, sample problems, and documentation (same as Appendix B) are 
available on a compact disk (CD). 

Variable- Curvature Shell Structures. A concerted effort has been expended to verify 
the formulation and implementation of the variable-radius shell analysis without complete 
success. Representation of the variable radius has been attempted using a segmented or 
superelement approach based on the tools developed previously under this grant. This 
approach required the development of an assembly procedure and special “joining” func- 
tions along segment junctures. Using this approach, the buckling response appears to 
be artificially “stiffened” perhaps because of the approximations used along the segment 
junctures. 

As an alternate approach, a global function representation of the shell radius has been 
attempted using a Legendre polynomial of order one. This would permit the modeling of a 
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shell with a linear change in curvature as a function of the circumferential coordinate. The 
present analysis method for buckling of anisotropic shells with variable curvature uses a 
segment approach where displacement fields within each segment are represented by Bezier 
polynomials and a first-order shear- deformation theory is used. In general, segments can 
be used in both axial and circumferential directions, however the present implementation 
considers only segments in the circumferential direction. Continuity of displacement at 
the junctures of adjacent segments are imposed using C° and C 1 conditions obtained from 
the properties of the Bezier control points. The shell with variable curvature is assumed 
to consist of two or more curved panels of constant curvature which is representative of 
fuselage or wing structures. 

Results are presented for a composite cylindrical panel subjected axial compression, 
a non-circular fuselage segment subjected to axial compression, and a composite wing 
leading-edge variable-curvature panel subjected to combined axial compression and shear. 
Sanders-Koiter shell theory is used in these studies. Buckling loads from the present anal- 
ysis are compared with those obtained from the STAGS finite element code. The STAGS 
finite element model consists of the 410 element, and curved surfaces are approximated 
as an assembly of flat surfaces. Buckling loads obtained from finite element solutions are 
determined to be four percent lower than those of the present analysis. Implementation 
of this method has been verified and documented (see Appendix C). These results were 
reported in a 1998 conference paper (CP9), and a journal paper will be submitted to the 
Composite Structures journal in the near future. 

The source code, sample problems, and documentation (same as Appendix C) are 
available on a compact disk (CD). 

Sandwich Plates and Shells. Sandwich construction techniques offer many advantages 
that can be exploited for advanced vehicles such as the HSCT. Recently PANDA2 has 
been extended to handle panels with sandwich wall construction by including the additional 
failure modes. Using this tool, various sandwich panel designs were assessed, and dominant 
designing failure modes based on the various mathematical models were identified. An 
exploratory study of the design of sandwich panels was performed using the PANDA2 
software system with the cooperation of a research scientist, Dr. David Bushnell. 

The status of analysis methods available was introduced, and analysis needs and 
refinements were defined. Assessment of sandwich panels, and their known potential failure 
modes and mechanisms will be performed using PANDA2 and their impact on the design 
process is identified. 

PANDA2 analysis is based on a global single layer approach wherein the sandwich core 
material is treated as just another layer in the laminate for determining global buckling 
behavior. Local analyses of failure modes account for core materials and the different 
face sheets, if applicable, in an analytical approach using solutions from Plantema, Vinson 
and PANDA2’s models. Failure modes include face wrinkling, face dimpling, core shear 
crimpling, core transverse shear stress failure, core crushing and tension, and face sheet 
pull-off. Of the approximations built into PANDA2, those associated with the transverse 
shear effects and the single-term buckling solution used in the PANDA-type (closed-form) 
analysis appear to be the more limiting factors in the analysis of sandwich structures. 

This work was accomplished by a Master’s graduate student (Mr. Hao Jiang) under 


9 



the direction of Dr. Knight (see T3). He should defend his thesis during the Summer 1998 
semester and a copy of his thesis will be forwarded to the technical monitor at that time. 
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APPENDIX A 


Documentation for the Analysis and Design Code 
for the Buckling of Grid-Stiffened Panels 


• Theoretical Foundations 

• User Instructions 

• Sample Problems 
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APPENDIX B 


Documentation for the Analysis and Design Code 
for the Buckling of Grid-Stiffened Cylinders 


• Theoretical Foundations 

• User Instructions 

• Sample Problems 
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Chapter 1 
INTRODUCTION 


The theoretical foundation and user instructions for a FORTRAN code for the de- 
sign and analysis of composite grid-stiffened circular cylindrical panels exhibiting global 
buckling are presented. Buckling analysis of composite grid-stiffened panel is performed 
by an analytical tool using an improved smeared stiffener theory ([1]) for the global 
buckling analysis, and a Rayleigh- Ritz-type buckling analysis for skin segment with gen- 
eral parallelogram-shaped ([2]) and general triangular planform ([3]) to assess local skin 
buckling. Crippling of stiffener segments are assessed by a method given in Reference [4]. 

The integration of this analysis method with a design optimization process for dis- 
crete design variable, such as the genetic algorithm ([5] is done to obtain a design opti- 
mization tool for grid-stiffened panel. The optimization tool ([6]) provides optimal de- 
signs for a buckling resistant grid-stiffened panel for a given set of in-plane design loads, 
boundary conditions of the panel, panel material properties, and the length, width, and 
radius of the panel. The design variables are the height and thickness of the stiffener, the 
axial and transverse stiffener spacings, the skin laminate and the stiffening configuration 
(isogrid, orthogrid, etc.). 

The theoretical foundation for the analyses involved in the buckling of grid-stiffened 
panel are discussed briefly in Chapter 2. Instructions for setting up input files for the 
FORTRAN code are given in Chapter 3. To provide flexibilities in performing different 
types of optimization, instructions are also provided in Chapter 3 for modifying the code 
to adjust the number of design variables to facilitate a particular type of optimization. 
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Chapter 2 

THEORETICAL FOUNDATIONS 


The theoretical foundations for the improved smeared stiffener theory ([1]), buckling 
of panels with general parallelogram and triangular-shaped planform ([2, 3]), and crip- 
pling of stiffener segments ([4]) are discussed in this chapter. The optimization strategy 
for buckling of grid-stiffened panel with global and local buckling constraints are also 
discussed in this chapter. 

2.1 IMPROVED SMEARED STIFFENER THEORY 

The improved smeared stiffener theory for stiffened panels, used here includes skin- 
stiffener interaction effects. Skin-stiffener interaction effects may lead to overestimation 
of buckling loads especially when the stiffener spacings are not small. 

If a stiffened plate is bent while it is supported on all four edges, the neutral surface 
in the neighborhood of the stiffener will lie between the mid-plane of the skin and the 
centroid of the stiffener. It is convenient to think of this as a shift of the neutral surface 
from the centroid of the stiffener. Hence, the approximate stiffness added by a stiffener to 
the skin stiffness will then be due to the skin-stiffener combination being bent about its 
neutral surface rather than due to the stiffener being bent about its own neutral surface 
or the skin neutral surface. The shifted location of the neutral surface is determined the- 
oretically through a study of the local stress distribution near the skin-stiffener interface 
for a panel with a blade stiffener. 

The neutral surface profile of the skin- stiffener combination is developed analytically 
using the minimum potential energy principle and statics conditions. The skin-stiffener 
interaction is accounted for by computing the bending and coupling stiffness due to the 
stiffener and the skin in the skin-stiffener region about a shifted neutral axis at the 
stiffener. 

A grid-stiffened panel may be considered to be an assembly of repetitive units or 
unit cells (see Figure 2.1). Any stiffener segment in the unit cell may be isolated in a 
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semi-infinite skin-stiffener model as shown in Figure 2.1 for a diagonal stiffener. The 
approach for obtaining the neutral surface in a semi-infinite stiffened panel is given in 
Reference [1]. 

A typical profile of the neutral surface for a skin-stiffener combination is shown in 
Figure 2.2. The distance y* represents the distance from the centerline of the stiffener 
to the point where the neutral surface coincides with the mid-surface of the skin. The 
average of the neutral profile over the distance y* is Z*. The quantities y* and Z* are 
obtained numerically. 

The smeared stiffnesses of a stiffened panel is obtained by mathematically converting 
the stiffened panel to an equivalent unstiffened panel (Ref. [7]. The smeared stiffnesses 
are developed on the basis that the strain energy of the stiffened panel should be the 
same as that of the equivalent unstiffened panel. These smeared stiffnesses can then be 
used in a Rayleigh- Ritz type analysis to solve for buckling loads of the stiffened panel. 
In Reference [7], the strain energy of the skin and stiffeners in the unit cell is obtained by 
using stiffnesses of the skin and the stiffeners which are computed about the mid-surface 
of the skin. Since, there is a shift in the neutral surface at the stiffener, the stiffness of 
the stiffeners and the skin segment directly above it has to be computed about a shifted 
neutral surface so as to account for the skin-stiffener interactions. 

The correction to the smeared stiffnesses due to the skin-stiffener interaction is herein 
introduced by computing the stiffness of the stiffener and the skin segment directly con- 
tiguous to it according to the following criteria. 

1. If y* < 2/4, then the reference surface for the stiffener is Z n . 

2. If y* > t/A, then the reference surface for the stiffener is Z*. 

In either case, the reference surface of the skin is taken to be its mid-surface. Other 
more elaborate and accurate schemes can be used to introduce the skin-stiffener interac- 
tion using the neutral surface profile. However, the one described herein is simple, and 
provides sufficiently accurate buckling loads for the preliminary structural design ([!]). 
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2.2 LOCAL BUCKLING OF SKIN SEGMENTS 

The shape of a skin segment on a grid-stiffened panel depends on the stiffening con- 
figuration. If the stiffening configuration involves diagonal stiffener only, then the skin 
segment has a rhombic planform. If the stiffening configuration has diagonal stiffen- 
ers with axial or transverse stiffeners, then the skin segment has an isosceles triangular 
planform. For a general grid-stiffened panel, the skin segment has a right-angle trian- 
gular planform, and for an isogrid panel the skin segment has an equilateral triangular 
planform. 

Buckling analyses for panel with these kinds of planforms is achieved through the 
use of ” circulation function” and accounts for material anisotropy, different boundary 
conditions, and combined in-plane loading. A First-Order Shear-Deformation Theory is 
used. The shell theory that is used can be either Sanders- Koiter, Love, or Donnell theory. 
This is achieved through tracer coefficients. 

2.2.1 Physical and Computational Domain 

The buckling analysis of these local skin segments is enhanced by mapping their 
physical domain into a computational domain. Consider a general quadrilateral or tri- 
angular panel subjected to a state of combined in-plane loading where the loading and 
material properties are defined using the coordinate system (x — y ) shown Figure 2.3. 
The transformation from a physical domain to computational domain is necessary when 
dealing with general quadrilateral and triangular geometries in order to facilitate the com- 
putation of linear stiffness and geometric stiffness matrices and imposition of boundary 
conditions. 

The physical domain V[x,y] is transformed to a computational domain Z>[£, 77 ] as 
indicated in Figure 2.3. The mapping for a quadrilateral is 

4 

i= 1 

y(Z,ri) = ( 1 ) 

2 = 1 

where X{(i = 1,2, 3,4) and yi(i = 1,2, 3,4) are the physical coordinates of the i th corner 
of the panel, £ and 77 are the natural coordinates for the quadrilateral geometries, and 
Ni ( i = 1 , 2 , 3, 4) are the bilinear mapping functions given by 
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Ni(t,v) = j(l -?)(! + V) 
m,v) = \( i+fKi+v) 
m,n) = \( 1+0(1 -I?) 

N 4 (t,y)= l -(\ - 0(1 - If) 
The Jacobian of the transformation is 


' dx 

dy - 

di 

d£ 

dx 

dy 

. dr] 

dr] . 


which is independent of the natural coordinates for general parallelogram-shaped geome- 
tries. This results in substantial computational savings in the overall formulation. 

The mapping for a general triangle is 

x(£,,Ti,p) = £xi -f rjx 2 + px z 

yttiViP) = £yi + vv2 + py3 (3) 

where (, p and p are the area coordinates for the case of triangular geometries, and 
Xi(i = 1,2,3) and yi(i = 1,2,3) are the physical coordinates of the i th corner of the 
panel. Note that the third area coordinate will be expressed in terms of the other two 
or p = (1 — f — 7]) based on the constraint that the sum of the area coordinates must be 
equal to one. The Jacobian of the transformation is independent of the area coordinates. 
The Jacobian, in either case, is used to relate derivatives in the two domains. 

2.2.2 The Rayleigh-Ritz Method 

The Rayleigh-Ritz method is an approximate method for solving a certain class of 
problems. Accordingly, trial functions with some unknown coefficients and satisfying the 
essential or geometric boundary conditions are introduced in the energy functional of the 
problem. The minimum conditions of this functional are then imposed, and resulting 
algebraic equations are solved for the unknown coefficients. These trial functions are 
called the “Ritz” functions. 

The Ritz functions used here are expressed in terms of natural coordinates for the 
quadrilateral geometry or area coordinates for the triangular geometry for displacement 
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field. The components of the displacement vector are three translations (D\, D 2 , D 3 = 
u 0 ,Vo,w) and two cross-sectional or bending rotations (D 4 ,D 5 = when consider- 

ing transverse-shear deformation effects. Each displacement component is approximated 
independently by a different Ritz function. The approximation for the i th component of 
the displacement vector is given by 

N 

= Yl a i3 d ij 

3 = 1 
N 

= for 2 = 1,2,3, 4,5 (4) 

3 = 1 

where dij represents the j th term in the TV- term approximation for the i th displacement 
component, ajj are unknown coefficients to be determined, and Ih(£, 77 ) are the circulation 
functions. 

The circulation functions T,- in Equation (4) are then used to impose different bound- 
ary conditions along each edge of the plate. Each term I\ is the product of three functions 
in the case of the triangular plate geometry and four functions in the case of the quadri- 
lateral plate geometry. Each function is the equation of an edge of the triangular or 
quadrilateral plate as shown in Figure 2.3 raised to an independent exponent for each 
displacement component. Thus, the circulation functions for the quadrilateral plate are 

r< = (i-i7)«(i-e') w (i+i/r(i+o j< 

and for the triangular plate are 


r, = w(i -€-«?)” (5) 

For example, considering the quadrilateral plate case, pi refers to edge 1 , qi refers to edge 
2 , r,- refers to edge 3, s,- refers to edge 4 as indicated in Figure 2.3. These exponents 
are used to impose different boundary conditions. If the i th displacement component is 
free on a given edge, then the exponent for that edge will have a value of zero. If the 
i th displacement component is constrained on a given edge, then the exponent for that 
edge will have a value of one. Only geometric boundary conditions are imposed in this 
approach. Thus, a simply supported condition for bending fields can be imposed on edge 
1 by setting: 

• p 3 = 1 for w, p 4 = 0 for <j> XJ p$ — 0 for (f> y 
A clamped condition for bending fields can be imposed on edge 1 by setting: 
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• P 3 = 1 for w, £>4 = 1 for <j> x , ps — 1 for <j) y 

A free-eclge condition can be imposed on edge 1 by setting: 

• pi — 0 for uq , Vo , w , 4> x and (j) y 

The term fj in Equation (4) is a polynomial function in £ and 77 , and in its simplest 
form is a power series in f and p and is expressed as 

fAU) = H m ’n n ’ 

m h n j = (0,0), (1,0), (0,1), (2,0), (1,1), (0,2),... (6) 

The values of mj and rij are used basically to define terms in a two-dimensional Pascal’s 
triangle. The number of terms N in Equation (4) defines the order of a complete function 
in two variables. The table below gives the value of N for polynomials of different degrees. 

Table 1 Degree of polynomials with value of N 


N 

Degree of 
polynomials 

N 

Degree of 
polynomials 

1 

0 

45 

8 

3 

1 

55 

9 

6 

2 

66 

10 

10 

3 

78 

11 

15 

4 

91 

12 

21 

5 

105 

13 

28 

6 

120 

14 

36 

7 

136 

15 


2.3 CRIPPLING OF STIFFENER SEGMENT 


The local stiffener segment is analyzed to determine whether stiffener crippling will 
occur. Reference [4] provides a method for determining the buckling load of a stiffener 
segment. Accordingly, the stiffener segment at the nodes or intersection points of stiff- 
eners are assumed to be clamped while the stiffener- skin attachment is assumed to be a 
simple support. From Ref. [4], the crippling load of the stiffener is N cr i p and is given by 


N ■ 

1 v crip 

il 

i i + y - ‘i 


N c i 

= z[ 

i^E n 


12I?[1 - (v? 2 E 22 /E u )} 


T 


£12 1 

h 2 J 


where 


( 7 ) 




where 


s z = |G r i 3 / 5 , is a shear correction factor, 

L\ = 2 L is the length of the stiffener, 
h is the height of the stiffener, 
t s is the thickness of the stiffener, 

En is the longitudinal modulus of the stiffener material, 

E 22 is the transverse modulus of the stiffener material. 

2.4 OPTIMIZATION OF GRID-STIFFENED PANEL 

The analysis and design of grid-stiffened composite panels subjected to combined 
loads require several key steps. In the present optimization procedure, acceptable designs 
are those which buckle globally and do not exhibit any local skin buckling or stiffener 
crippling. The first step is to assess the global buckling response of the grid-stiffened 
panel. Once this global buckling response is determined, the second step is to determine 
the local skin buckling response for general quadrilateral and/or triangular skin segments 
that occur locally between stiffeners. The third step is to determine whether stiffener 
buckling or stiffener crippling has occurred at this global load level. This sequence of 
steps is performed repeatedly in a design cycle until an optimum or near-optimum design 
is obtained. The genetic algorithm is used herein and the buckling analyses involved in 
the global and local buckling of grid-stiffened panels have been discussed in Sections 2.1 
and 2.2. 


2.4.1 Panel Design Procedure 

The design of grid-stiffened composite panels requires that many of the design vari- 
ables, such as stiffener spacing and stiffener thicknesses may only take on certain discrete 
values rather than vary continuously over the design space, and often a “family” of good 
designs is needed rather than a single-point design due to manufacturing requirements. 
Gradient-based methods for structural optimization are not appropriate in this case. 

The genetic algorithm is a method for “evolving” a given design problem to a family 
of near-optimum designs ([5]). Based on Darwin’s theory of survi val-of-the- fittest, the 
genetic algorithm involves the random creation of a design population that “evolves” 
towards some definition of fitness. The genetic algorithm is attractive due to their sim- 
plicity of approach in discrete variable combinatorics. The genetic algorithm can be used 
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directly to solve unconstrained optimization problems, while constrained optimization 
must first be transformed to an unconstrained optimization problem (e.g., use of an ex- 
terior penalty function). Stochastic processes are used to generate an initial population 
of individual designs and the process then applies principles of natural selection and 
survival of the fittest to find improved designs. Furthermore, since the discrete design 
procedure works with a population of designs it can explore a large area of the design 
space and climb different hills. This is a major advantage as the converged solution 
contains many optima of comparable performance. The cost of having a large number of 
function evaluations is offset by the fact that a large number of optima solutions are now 
available. In a gradient-based optimization procedure, only a single-point design, usually 
the extremum to the starting point, is obtained. The genetic algorithm produces a pop- 
ulation or family of good designs which may include the global optimal design, rather 
than a single design. Hence, it is an appropriate tool for designing general grid-stiffened 
panels. 

2.4.2 Design Problem Definiton 

The present design problem is to minimize the weight per unit area of a grid-stiffened 
composite panel given the design loading condition, the length and width of the panel, 
the material properties for the skin and stiffeners, and the boundary conditions of the 
panel. The design variables include stiffener spacings (a, b ), the stacking sequence of 
the skin, stiffener layout or stiffening configuration, stiffener thickness (i a ), and stiffener 
height (hi = h 2 — h 3 = h) as shown in Figure 2.4. The axial and transverse directions of 
the panel are along the x and y axis respectively. All stiffeners are assumed to be of the 
same height and thickness for manufacturing and assembly reasons. The design sought 
here is a panel of minimum weight in a certain design space which buckles globally at 
the design loads. The design is defined by setting up the optimization procedures in the 
following way. First, the global buckling load is assumed to be a scalar multiple of design 
loads and has the form 


N x — Ao-Ah, N y — A gAT 2 , N xy — \qN\ 2 (8) 

where Ah, N 2 , N\ 2 are the applied in-plane prebuckling loads. These values represent 
the design loads for the grid-stiffened panel. Second, the design constraints imposed on 
panel include 
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1. The critical buckling load should be greater than or equal to the design loads, that 
is, A a > 1. 


2. Skin segments should not buckle at the critical buckling load, that is, \ s k > 1. 

3. Stiffener segments should not cripple at the critical buckling load, i.e., A 1} A 2 , A 3 > 
1 where Ai, A 2 , A 3 is the crippling load factor of the x-direction stiffener, y- 
direction stiffener and diagonal stiffener, respectively. 


The general form of each constraint equation is written as 


Pi = (t — i) < °-o i = 1, N c 

* 3 

Finally, the “Fitness” expression based on exterior penalty function approach is 


Fitness = ( 


Q 


F(X,n) 


) = Max- 


Q 


VK(X) + r,Ef [|ft(X)|+ a (X)]» 


( 9 ) 


( 10 ) 


where 

X = design variable vector 

i^X, rj) — modified objective function 

VK(X) = weight of panel per unit area 

r % Ef c [bi(X)| + ^(X)] 2 = penalty function 

Q = normalizing constant 

N c = number of design constraints 

r; = penalty parameter 

i — generation or iteration cycle in the optimization procedure. 

Once the global buckling load factor has been determined using the improved smeared 
stiffener theory, the loads acting on the stiffener and skin segments have to be deter- 
mined by distributing the loads based on the extensional stiffness of the skin and the 
stiffener. The procedure for distributing the applied loads for a general grid-stiffened 
panel involves three steps. First, the extensional stiffness coefficients for grid-stiffened 
panel are computed as follows (Ref. [7]): 


(-Aii)t 


h 2(Ah) 3 h sin 3 $ 

i + l 

2(An) 2 h t 2(Ah) 3 h cos 3 9 
r 

a a 

2 (An )3 h cosO sin 2 9 

f- (yi 66 ) 


+ (^ll)s 

+ (^22)5 


(^22)7 = 
{Aqq)t = 


a 


( 11 ) 
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where (Ah)t is total smeared axial extensional stiffness of the grid-stiffened panel, ( A^t 
is the total smeared transverse extensional stiffness of the grid-stiffened panel, (A 6 6 )t is 
the total smeared in-plane shear stiffness of the grid stiffened panel, (dn)i, (^ 11 ) 2 , (^ 11)3 
are the extensional stiffness of the axial, transverse and diagonal stiffeners, respectively, 
( Aij) s is the extensional stiffness of the skin, 0 is the orientation of the diagonal stiffener, 
and h is the height of the stiffener. Second, the loads carried by the skin segment which 
could be either a general parallelogram-shaped geometry or a general triangular-shaped 
geometry, at the panel global buckling load are 


(Nx)sk = 
( Ny )sk ~ 
(Ahy)sfc — 


(^11)3 

{A\\)t 

(^ 22)3 

(^22)t 

(^.66)3 

{Aqq)t 


N x = 

Ny = 

N xy = 


(^ll)s 

(Ah)t 

(^-22)3 

(^22)t 

(^d-66). 


(^•66)' 


A G iVi 
A G A ^ 2 
L A G A^12 


( 12 ) 


These values then correspond to the design loads used for the in-plane prebuckling load 
in the skin-segment local buckling computation. If the critical buckling load factor of the 
skin segment X s k is greater than or equal to one, then the skin-segment buckling load is 
greater than or equal to the global buckling load of the grid-stiffened panel. Third, the 
loads carried by each stiffener are computed. The load carried by the axial stiffener is 


(Ay 1 = 


(^n)i 


A G Ah — X\N C rip 


(13) 


where N cr { p is determined using Equation (7). The critical buckling load factor, Ai, of the 
axial stiffener has to be greater than or equal to one. The load carried by the transverse 
stiffener is 

(N x h = EiT.A G iV 2 = A 2 N cr , p (14) 

\A22)T 

and the critical buckling load factor, A 2 , of the transverse stiffener has to be greater 
than or equal to one. The load in the diagonal stiffeners has components from the axial, 
transverse, and in-plane shear loadings and is given by 


i^N — * Nd x sm0 T -f- (A ^dxy'jxCOsQ T (ATr^ySZTi^ — XsNcrip 


A fdx 


Ndy 


(Au)3sin 3 9 

(^ii)t 

(A n ) 3 cos 3 0 

(A<22)t 


A g Ai 


A G A r 2 


where 
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( Ndxy )r 
(Ndxy)y 


( A\\)2,cos9sin 2 0 
(A Q q)t 

(A u ) 3 cos0sin 2 0 
( A 66 )t 


-^ gN 12 
a 


XgNw 


( 15 ) 


Ndx is the contribution from the axial in-plane loading, Nd y is the contribution from the 
transverse in-plane loading, (Nd xy ) x is the contribution from the in-plane shear loading 
along the edge where x is constant, and (Nd xy ) y is the contribution from the in-plane 
shear loading along the edge where y is constant. The critical buckling load factor, A 3 , 
of the diagonal stiffener has to be greater than or equal to one. 


The weight per unit area of the grid-stiffened panel is 


W 

where 

Wi 

w 2 

w 3 

w s 


= A[(wi +W2+W3 + W s ) 
ab 

— 2 h a t s 
= 2 h b t s 

= 2 h t y/ a 2 + 6 2 

— d b tgfcn 


(16) 


w\ is the volume of the axial stiffeners in the unit cell, w 2 is the volume of the transverse 
stiffeners in the unit cell, w 3 is the volume of the diagonal stiffeners in the unit cell, w s 
is the volume of the skin in the unit cell, t s ki n is the thickness of skin, and p is the mass 
density of the material. 

2.4.3 Design Process Based on Genetic Algorithm 


Implementation of the genetic algorithm is shown schematically in Figure 2.5. The 
design process begins with a random selection of a specified number of designs which 
comprise the initial population (i.e., first generation) for the genetic algorithm. Material 
properties, length and width of panel, boundary conditions of the stiffened panel, and 
design loadings are input to the analysis processor routine. The buckling analysis is 
performed which provides the critical eigenvalues for the global buckling response of the 
grid-stiffened panel and the local buckling response of the skin and stiffener segments, 
which also computes the weight per unit area of the grid-stiffened panel. This procedure 
is repeated for each design configuration in the population. The “fitness” processor then 
evaluates the “fitness” of each design using Equation (10) and assigns a rank based on the 
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fitness expression or objective function. The current population of design configurations is 
then processed by the genetic operators (crossover, mutation, and permutation) to create 
a new population of design configurations for the next generations which combine the 
most desirable characteristics of previous generations. Designs from previous generations 
may be replaced by new ones (i.e., children) except for the “most fit” designs (i.e., 
parents) which are always included in the next generation. The process is repeated until 
design convergence is obtained, which is defined herein by specifying a maximum number 
of generations ( NSTOP ) that may occur without improvement in the best design. 



Figure 2.1 Semi-infinite plate model for skin-stiffener element. 



Figure 2.2 Typical profile for skin-stiffener element neutral surface. 
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Figure 2.3 Physical and computational domain for plate geometries. 
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(D Diagonal Stiffener 

Figure 2.4 Unit cell showing design variables. 









Figure 2.5 Flow chart showing the design optimization process. 
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Chapter 3 

USER INSTRUCTIONS 


User instructions for using two FORTRAN codes are provided in this chapter. The 
first code is for analyzing a grid-stiffened panel subjected to combined in-plane loading. 
The code provides the global buckling load, and local buckling load of skin and stiffener 
segments as output. Instructions for using this code are given in Section 3.1. The second 
code is for optimizing a grid-stiffened panel design subjected to global and local buckling 
constraints. Instructions for using this code and modifying it to obtain a particular type 
of optimization are given in Section 3.2. 

3.1 ANALYSIS CODE 

The analysis code is found in directory “panel/grid” and a makefile is used to link all 
the subroutines together. The makefile may have to be modified to account for different 
computer system. Grid-stiffened panel with the unit cell geometry as shown in Figure 
2.4 can be analyzed by using the code. Skin segments of the grid-stiffened panel are 
assumed to be simply-supported. Other boundary conditions of the skin segment can be 
accommodated through simple modifications to the source code. The executable for this 
code is “run” as specified in the makefile. 

3.1.1 Examples for Input and Output file 

An input file for a flat grid-stiffened panel with axial and diagonal stiffeners is given. 
The panel is 20.0 in. long and 56.0 in wide, and has an axial and transverse stiffener 
spacings of 3.3333 in. and 5.89473 in., respectively. The height and thickness of the 
stiffener is 0.5 in. and 0.06 in., respectively. The skin laminate has a ply stacking 
sequence of [60/0/60] with ply thickness of 0.006 in. The ply orientations are measured 
from the x-axis in the counter-clockwise directions. The stiffener is made of unidirectional 
material. The material for the skin and stiffener is assumed to have the following nominal 
ply mechanical properties; E\\ = 20.2 Msi; E 22 = 1.9 Msi; G \2 = G 13 = G 23 = 0.73 Msi 
and V \2 = 0.3. The panel is simply supported and is subjected to axial compression 
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loading of N x ~ 400 Ibs/in. The input file is named “pan.inp” and is given in List 1. 
Text after the character “!” are comments and need not be included in the actual file. 

Some considerations for the input file are listed below. 

1. The maximum number of plies in a laminate is 50, and the maximum number of 
material is 5. 

2. When a stiffener type (axial, transverse, or diagonal) is not present, its thickness 
and ply thickness are entered as zero. But not the material properties and, its 
height as shown for the transverse stiffener in the above input file. 

3. The dimension of the panel is read as (x, y) coordinate of each node (see Figure 
2.3). The radius of the panel is read by specifying the radius at each node. For a 
flat panel, the radius is input as a very large number. 

4. When analyzing panels stiffened in one direction only, one of the stiffener spacing 
is redundant. For example, an axially stiffened panel will have its stiffener spacing 
specified by the width of the unit cell only. The length of the unit cell is entered 
as the length of the panel. 

5. The number of terms or modes (N) used in the analysis is taken from Table 1. 
The maximum value for N is 150 as determined by parameter “nmod” in the file 
“panel. inc”. Using a value of N between 55 and 78 is usually sufficient. 

Finally, if the user wishes to change the boundary conditions of the skin segment in the 
analysis, the subroutine “bclocal.f” has to be modified. The output file is “pan. out” and 
is given in List 2. 

3.2 OPTIMIZATION CODE 

The optimization code will optimize a grid-stiffened panel for minimum weight sub- 
jected to global and local buckling constraints, and is found in directory “panel/optimize”. 
The executable for this code is “run” as specified in the makefile. The design variables 
are 


1. Axial stiffener spacing (a). 

2. Transverse stiffener spacing (6). 
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3. Stiffener height (h). 

4. Stiffener thickness ( t s ). 

5. Skin laminate (LAMI). 

6. Stiffening configuration (IGEO). 

The number of design variable is defined by the parameter “n” in the main program 
“main.f’ Each design variable can assume eight discrete values as allowed by the FOR- 
TRAN code. The eight discrete values of each design variable define the design space 
for optimization. The discrete values for os, 6, A, and t s are supplied through the file 
“inp.gen” which is read by the main program “main.f”. Part of the main program where 
the parameter “n”, and the parameter for the population size “m” are defined, the values 
for a, b , h, and t s are read is shown in List 3. The discrete values for LAMI and IGEO 
are given in Table 2. 


Table 2 Design space for design variables ICON and LAMI. 


Integer 

LAMI 

IGEO 

value 



1 

[±45/0] 2s 

axial stiffeners 

2 

[±45/90j 2s 

transverse stiffeners 

3 

[±45/0/90] 23 

axial and transverse stiffeners 

4 

(±45/0 2 ] 2a 

diagonal stiffeners 

5 

[±45/90 2 ] 2s 

axial and diagonal stiffeners 

6 

[±45/0 2 /90] 2s 

transverse and diagonal stiffeners 

7 

[±45/0/90 2 ] 2s 

axial, transverse and diagonal stiffeners 

8 

[±45/0 2 /90 2 ] 2j 

no stiffeners 


The laminate stacking sequence corresponding to various discrete values of LAMI 
are hard-wired in subroutine “cskin.f’, and can be changed by modifying subroutine 
“cskin.f’. The ply thickness in subroutine “cskin.f’ is kept constant for all laminates, and 
is read in “cskin.f’. Subroutine “cskin.f’ can be modified by the user to accommodate 
various laminate stacking sequences. The discrete values for IGEO are assigned in 
subroutine “panel.f’ and part of the code where IGEO are being assigned is shown in 
List 4. Subroutine “geom.f’ assigns the stiffening configuration based on the value of 
IGEO which is supplied by the main program “main.f’. 

Some parameters that may affect the optimization process are 
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• The population size “m” is hard-wired in “main.f’. Usually m > 2n, and this 
condition has been found to work well, and m 2n is not recommended. 

• The probabilities of crossover, mutation, and permutation have been hard-wired 
to 1.0, 0.1, and 0.95 in “main.f’ . These values work very well for the optimization 
problem under consideration. 

• The termination criteria “NSTOP” is hard- wired in “main.f’. The user has to 
experiment with the value of “NSTOP”. Usually the code is run with a value of 
“NSTOP” and then with another value of “NSTOP” greater than the previous one. 
If there is no change in the optimal designs, then the second value of “NSTOP” 
provide a good value as a stopping criteria. For the problem under consideration, 
NSTOP=25, is usually sufficient. 

• The penalty parameter r; in Equation (10) can either increase at every i th gen- 
eration or can be constant for all generation. In subroutine “panel. f’, r t - is kept 
constant at 1000. By commenting the line where “ ainipen = 1000.0”, the user 
can set 

r,- = 1000 + i 2 (17) 

Keeping r* constant works very well for the present optimization problem. 

According to Lists 3 and 4, the code has been set up to optimize grid-stiffened panel 
with all design variables active. 

3.2.1 Changing the Type of Optmization 

The user may want to optimize a grid-stiffened panel with less number of design 
variables. For example, the skin laminate, and the stiffening configuration may be fixed, 
and the only design variables are a, 6, /i, and t s . An example of such an optimization is 
provided with the required input files, and the output files from the code are explained. 

Consider the panel described in Sub-section 3.1.1, the panel is to be optimized for N x 
— 400 lbs/in., with design variables being a, 6, h, and t 3 . The skin laminate is [60/0/ 60]^ 
with a ply thickness of 0.006 in. Only axial and diagonal stiffeners are considered and 
therefore I GEO = 5. The axial stiffener spacing a and transverse stiffener spacing b is 
treated as one design variable i.e., (a, b) is a design variable. Values for (a, b) are provided 
such that the stiffening configuration closely approximates an isogrid configuration (i.e., 
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9 30°). Hence, there are three design variables. List 5 and 6 show the appropriate 

modifications to “main. f 5 and ” panel. f” respectively. In List 5, “n” has been changed to 
3, and “m” has been changed to 8. While in List 6, a “! modify” indicates the line that 
has been modified. 

The code needs two input files, namely “inp.gen”, and “pan.inp”. The file “inp.gen” 
is read by program “main.f” and it defines the design space for the stiffener spacings, 
and the height and thickness of the stiffener. The file “pan.inp” provides the problem 
parameters for the optimization problem and is read by subroutine “panel. f’. Example 
for “inp.gen”, and “pan.inp” are given in List 7 and 8 respectively. 

The output files produced by the code are “best. gen”, “on. gen”, and “pan. out”. The 
files “best. gen”, and “on. gen” are produced by program “main.f’, and the file “pan. out” 
is produced by subroutine “panel. f”. The optimal designs ranked according to Equation 
(10) are stored in the file “best.f” and the convergence history of the optimization is 
stored in the file “on. gen”. The file “best. gen” and “on. gen” for the above example is 
given in Lists 9 and 10 respectively. 

The file “pan. out” stores the information about each design resulting from the anal- 
ysis and is quite large. It is useful in obtaining the buckling loads and other information 
about the optimal designs stored in file “best .gen”. To access information about an opti- 
mal design, use the value of its fitness (FS) in “best. gen” and locate that number (critlb) 
in “pan. out” using the search option of the unix editor being used. For example, the best 
design, which is the first design in “best. gen” has “FS= 0.7669355E+03” . Searching for 
the pattern “0.7669355E-H03” in “pan. out” will bring the cursor to where information 
about the best design is written. List 11 and 12 give information about the first and 
second optimal designs which have been extracted from “pan. out”. 
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0.036 ! 

1 ! 

1,20. 2e6 , 1 . 9e6 ,0 
6 ! 

1,1,0.006,60.0 ! 
2,1,0.006,0.0 

3.1.0. 006,-60.0 

4.1.0. 006,-60.0 

5.1.0. 006.0.0 

6.1.0. 006.60.0 

0.060 ! 

1 ! 

1,20. 2e6 , 1 . 9e6 ,0 
1 ! 

1,1,0.060,0.0 ! 

0.0 ! 

1 ! 

1,20. 2e6 , 1 . 9e6 , 0 
1 ! 

1 , 1 , 0 . 0 , 0.0 ! 

0.060 ! 

1 ! 

1,20. 2e6 , 1 . 9e6 ,0 
1 ! 

1,1,0.060,0.0 ! 

0.5, 0.5, 0.5 ! 

0 . 0 , 0 . 0 , 1 . 0el2 

20.0. 0.0.1. 0el2 

20.0. 56.0.1. 0el2 
0.0,56.0,1. 0el2 

1 . 0 . 1.0 
0 , 1 , 0,1 
1 , 1 , 1,1 
0 , 0 , 0,0 
0 , 0 , 0,0 

3.333333,5.89473 

55 

400.0,0.0,0.0 


LIST OF FILES 


List 1: Example for Input file (pan.inp) 


thickness of skin 
Number of material 

.3,0.73e6 ! Material No., Ell, E22, vl2, G12 
No. of plies 

layer No., Material No., ply thickness, theta 


axial stiffener thickness 
Number of material 

.3,0.73e6 ! Material No., Ell, E22, vl2, G12 
No. of plies 

layer No., Material No., ply thickness, theta 
transverse stiffener thickness 
Number of material 
.3,0 .73e6 ! Material No. 

No. of plies 
layer No., Material No. 
diagonal stiffener 
Number of material 
. 3,0.73e6 ! Material No. 

No. of plies 
layer No., Material No. 
height of axial, transverse, and diagonal stiffener 


Ell, E22 , vl2 , G12 
ply thickness, theta 


Ell, E22 , vl2 , G12 
ply thickness, theta 


(x,y), radius for node 1 for 20 in 

(x,y), radius for node 2 for 20 in 

(x,y) , radius for node 3 for 20 in 

(x,y) , radius for node 4 for 20 in 

B.C's of U for panel 
B.C's of V for panel 
B.C's of W for panel 
B.C's of Px for panel 
B.C's of Py for panel 
length, width of unit cell 
# of modes considered (N) 

Nx, Ny, Nxy (loading condition) 


by 56 in 
by 56 in 
by 56 


in 


by 56 in 


panel 

panel 

panel 

panel 
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List 2: Example of output file (pan. out) 


SKIN LAMINATE DATA 

STACK THICKNESS= 3 . 5999999999999997E-02in 


NO. of 

MATERIAL types^ 

: 1 



MAT. NO 

El 

E2 

V12 

G12 


(psi) 

(psi) 


(psi) 

1 

0 . 2020E+08 

0 . 1900E+07 

0.3000 

0 . 7300E+06 

LAYER 

MATERIAL 

THK 

ORIENTATION 


No. 

No. 

(in) 

(deg) 


1 

1 

0.0060 

60.0000 


2 

1 

0.0060 

0.0000 


3 

1 

0.0060 

-60.0000 


4 

1 

0.0060 

-60.0000 


5 

1 

0.0060 

0.0000 


6 

1 

0.0060 

60.0000 



X-STIFFENER 

LAMINATE DATA 


STACK THICKNESS 3 5 . 9999999999999998E-02in 


NO. of 

MATERIAL types= 

* 1 



MAT. NO 

El 

E2 

V12 

G12 


(psi) 

(psi) 


(psi) 

1 

0 . 2020E+08 

0 . 1900E+07 

0.3000 

0 . 7300E+06 

LAYER 

MATERIAL 

THK 

ORIENTATION 


No. 

No. 

(in) 

(deg) 


1 

1 

0.0600 

0.0000 



Y-STIFFENER LAMINATE DATA 

STACK THICKNESS- 0 . OOOOOOOOOOOOOOOOE+OOin 
NO . of MATERIAL types= 1 


MAT. NO. 


LAYER 

No. 

1 


El 
(psi) 

0.2020E+08 

MATERIAL 

No. 

1 


E2 

(psi) 


V12 


G12 

(psi) 


0 . 1900E+07 0.3000 


0 . 7300E+06 


THK 

(in) 

0.0000 


ORIENTATION 

(deg) 

0.0000 


DIAGONAL-STIFFENER LAMINATE DATA 

STACK THICKNESS 3 5 . 9999999999999998E-02in 
NO. of MATERIAL types 3 1 

MAT. NO. El E2 V12 G12 
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(psi) (psi) (psi) 


1 

0 . 2080E+08 

0 . 1880E+07 

0.3000 0 . 7400E+06 

LAYER 

MATERIAL 

THK 

ORIENTATION 

No. 

No. 

(in) 

(deg) 

1 

1 

0.0600 

0.0000 


Height of X-stiffener = 0.5000000000000000 
Height of Y-stiffener = 0.5000000000000000 
Height of Dia-stiff ener = 0.5000000000000000 


Node # X 


Y RADIUS 


1 

2 

3 

4 


0.000 

20.000 

20.000 

0.000 


0 . 000 1000000000000 . 000 
0.000 1000000000000.000 
56.000 1000000000000.000 

56.000 1000000000000.000 




BOUNDARY 

CONDITIONS 

Edge 

U 

V 

W 

Px 

1 

1 

0 

1 

0 

2 

0 

1 

1 

0 

3 

1 

0 

1 

0 

4 

0 

1 

1 

0 


Py 


0 

0 

0 

0 


STIFFENER ORIENTATION (deg) 
UNIT CELL LENGTH (in) 

UNIT CELL WIDTH (in) 


29.48715171286544 

3.333333000000000 

5.894730000000000 


No of MODES CONSIDERED = 


55 


* U * * V * * 

0 0 0 0 0 

10 10 1 
0 10 10 
2 0 2 0 2 

11111 


0 9 0 9 0 

LOADING MATRIX 

400.000 

0.000 


W * *pX * *pY * 

0 0 0 0 0 

0 10 10 
10 10 1 

0 2 0 2 0 

11111 


9 0 9 0 9 


0.000 

0.000 


END OF INPUT DATA 


SKIN STIFFNESS DATA 

Extensional Stiffness (lbs) 

319210.997 102680.460 0.000 

102680.460 319210.997 0.000 



0.000 


108265.268 


0.000 


Coupling Stiffness (lbs) 


0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

Bending Stiffness 

(lbs-in) 


29.504 

12.073 

5.245 

12.073 

37.478 

15.469 

5.245 

15.469 

12.676 

Transverse shear 

stiffness (lbs/in) 


0 . 21900E+05 

0 . 00000E+00 


0 . 00000E+00 

0 . 21900E+05 



STIFFENER EXTENSIQNAL COUPLING BENDING 
STIFFNESS (lbs) (lbs in) (lbs in~2) 


1 0 . 6060E+06 

2 0 . 0000E+00 

3 0 . 6240E+06 
« ONLY AXIAL & DIAGONAL 
Stiffening Parameter (X) 
Stiffening Parameter (Y) 
Stiffening Parameter (D) 


- . 1624E+06 

0 . 5615E+05 

0 

0. 0000E+00 

0. 0000E+00 

0 

- . 1672E+06 

0 . 5781E+05 

0 

STIFFENERS 

» 



= 0.7183978195291273 

= O.OOOOOOOOOOOOOOOOE+OO 
= 0.2029944179752305 


X-stiff ener 

znn = -0.1288619996825681 

zstar = -1 . 3159971231018306E-02 

ystar = 1.616924400000000 

nstep = 56 


D-stiff ener 

znn = -0.1294158431233926 

zstar = - 1 . 504842384267 1259E-02 

ystar = 3.292702233808636 

nstep - 99 


CORRECTED STIFFNESS 


STIFFENER EXTENSIQNAL COUPLING BENDING 
STIFFNESS (lbs) (lbs in) (lbs in~2) 


1 

2 

3 


0 . 6060E+06 

- . 1547E+06 

0 . 5198E+05 

0 

0. 0000E+00 

0. 0000E+00 

0. 0000E+00 

0 

0 . 6240E+06 

- . 1581E+06 

0 . 5292E+05 

0 


EXTENSIQNAL SMEARED STIFFNESS 


Stiffeners 

230840.070 

78957.215 

0.000 

stiffeners + skin 
550051.067 


78957.215 

246923.403 

0.000 

181637.676 


0.000 

0.000 

78957.215 

0.000 


SHEAR 
(lbs) 

. 2190E+05 
. 0000E+00 
2220E+05 


SHEAR 

(lbs) 

2190E+05 

OOOOE+OO 

2220E+05 
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181637.676 566134.400 0.000 

0.000 0.000 187222.484 

COUPLING SMEARED STIFFNESS 

Stiffeners 

-58876 .790 
-20008.824 
0.000 

stiffeners + skin 
-58876.790 
-20008.824 
0.000 

BENDING SMEARED STIFFNESS 
Stiffeners 

19776.506 6696.980 0.000 

6696.980 20943.508 0.000 

0.000 0.000 6696.980 

stiffeners + skin 

19806.010 6709.053 5.245 

6709.053 20980.986 15.469 

5.245 15.469 6709.656 

SMEARED TRANSVERSE SHEAR STIFFNESS (lbs/in) 

stiffeners 

11137.905 
0.000 

stiffener + skin 

11137.905 
0.000 


BEGIN BUCKLING ANALYSIS 


Global Lamda = 1.007544806578725 

(Nx) _sk = 233.0447616099682 

(Ny) _sk = O.OOOOOOOOOOOOOOOOE+OO 
(Nxy)_sk = 0.0000000000000000E+00 

(Nx).l = 892.3931403986120 

Ndx_l = 109.5527510697837 

Ndx_2 = 0.0000000000000000E+00 

(Ndxy) _x = O.OOOOOOOOOOOOOOOOE+OO 
(Ndxy) _y = O.OOOOOOOOOOOOOOOOE+OO 
DETERMINANT = 9.824549017545001 

skilara = 1.024367601996003 

riblx = 2.062871407188981 

ribld = 34.25216887196387 

Volume (lbs/ft ~2) (rho ~ 0.057 lbs/in~3) = 0.5487635861376152 

CPU TIME 3.071758 


0.000 

11594.610 

0.000 

11594.610 


-20008.824 

-62573.722 

0.000 

-20008.824 

-62573.722 

0.000 


0.000 

0.000 

-20008.824 

0.000 

0.000 

-20008.824 


o o o o o 
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List 3: Part of main program “main.f 1 


C GENETIC ALGORITHM 

C 

C IS(I,J), I being the individual number and J its Jth bits. 

C IS(I,J) bit string number I (old generation) 

C JS(I,J) bit string number I (new generation) 

C CRITLB(I) fitness associated to the individual I 

C FNS(I) normalized fitness of the individual I 

C M population size 

C NLA maximum number of layers 

C N number of bits in a string (=NLA/4) 

C PC probability of implementing crossover 

C PM probability of implementing mutation 

C PP probability of permutation 

C PRI probability of inversion 

C ITER iteration (generation) number 

C 
C 

parameter (n=6 ,m=14,nn=m*2) 
c 

c n = number of design variables 

c m = number of design in each group 

c nn= just for dimension ( parameter) 
c 

DIMENSION IS(nn,n) , FNS(IOO), RD(n) ,JS(m,n), 

&CRITLB(100) ,CRITZ(100) ,aas(8) ,bbs(8) ,hhl(8) ,tthkl(8) , 

&ISK(50 ,m) ,FSK(50) ,NGEN(50) ,alength(n) 
real*4 tcp2(2) 
c common/genpara/n,m,nn 

c common /stap/alength,area 

COMMON/PIE/PI 
COMMON/MATGEO/T , NLA 
common /function/critlb 
c 

OPEN (UNIT=15 , FILE='on.gen') 

OPEN (UNIT=12, FILE='best.genO 
OPEN (UNIT=9 , FILE- * inp .gen' ) 

C 

c nla=16 

c N=NLA/4 

n=15 


Maximum number of generations 

read (9,*) (aas(ij) ,ij=l,8) 
read (9,*) (bbs(ij) , ij=l,8) 
read (9,*) (hhl(ij) ,ij=l,8) 
read (9 ,*) (tthkl(ij) , i j = 1 ,8) 

C 

LTT=300 

C 



noo o o o o oooon oooon o o 
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Genetic parameters 

PC=1 . 00D+00 
PM=0 . 10D+00 
PP=0. 95+00 
M=16 


NFT is the number of evaluations of the objective function 
without improvement before the search stops. 

NFT=150 


Initialization of the stopping criterion 

NCRIT=0 
0PTI=0 . D+00 

NSTOP is the maximum number of generations without 
improvement . 

NST0P=15 

Initialization of the parameters of the subroutine STORE 
before the first call. 


call dtime(tcp2) 

write(12 ,*) 'CPU TIME =\tcp2(l) 

CL0SE(12) 

END 
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List 4: Part of subroutine “panel. f” 

subroutine panel (io , is , critlb , ainipen,nn,nd,aas ,bbs ,hhl , 
& tthkl) 


include "opt.inc" 
include "panel. inc" 
dimension is(nn,nd) 
c 

c ---------------- 

c UNIT 5 FOR THE INPUT DATA FILE 

c ---------------- 

open (5, f ile='pan.inp' ) 
rewind 5 
c 

c ------------------ 

c UNIT 6 IS FOR THE OUTPUT DATA FILE 

open ( 6 , f ile= ' pan . out ' ) 


c 

c 

c 

c 

c 

c 

c 

c 


c 


clen = aas(is(io,l)) 
cwid = bbs(is(io,2) ) 
hi = hhl(is(io,3) ) 
h2 = hi 
h3 = hi 

thick = tthkl(is(io,4)) 
lami = is(io,5) 
igeo = is(io,6) 

call geom(sthl , sth2 , sth3 , igeo , thick) 

write(6 , *) ’ 

write (6,*) 'Laminate = ; ,lami 
write(6,*) 'Stiff ener thickness = ', thick 
write (6,*) 'Stiff ener height = ' ,hl 


[1] READING ALL INPUT DATA FOR A LAMINATE 
USING SUBROUTINE ISKIN & STSKIN 


call i skin ( sthk , nmsk , exsk , eysk , vsk , gsk , amat sk , nsk , 
& lnosk,msk,tsk,thesk) 

call cskin(sthk,nmsk, exsk, eysk, vsk, gsk, amat sk, nsk, 
& lnosk,msk,tsk,thesk,lami) 


write(6,51)critlb(io) 
format ( 'critlb = 'e!4.7) 

write(6,*) ' 

return 

end 


51 



o o o o o 
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List 5: Modifications to main program “main.f” 

C GENETIC ALGORITHM 

C 

C IS(I,J), I being the individual number and J its Jth bits. 

C IS(I,J) bit string number I (old generation) 

C JS(I,J) bit string number I (new generation) 

C CRITLB(I) fitness associated to the individual I 

C FNS(I) normalized fitness of the individual I 

C M population size 

C NLA maximum number of layers 

C N number of bits in a string (=NLA/4) 

C PC probability of implementing crossover 

C PM probability of implementing mutation 

C PP probability of permutation 

C PRI probability of inversion 

C ITER iteration (generation) number 

C 
C 

parameter (n=3 ,m=8 ,nn=m*2) 
c 

c n = number of design variables 

c m = number of design in each group 

c nn 55 just for dimension ( parameter) 
c 

DIMENSION IS(nn,n), FNS(IOO), RD(n) ,JS(m,n), 

&CRITLB(100) ,CRITZ(100) ,aas(8) ,bbs(8) ,hhl(8) ,tthkl(8) , 

&ISK(50,m) ,FSK(50) ,NGEN(50) , alength(n) 
real*4 tcp2(2) 
c common/genpara/n,m,nn 

c common /stap/alength, area 

COMMON/PIE/PI 
COMMON/MATGEO/T , NLA 
common /function/critlb 
c 

OPEN (UNIT=15, FILE= i on.gen > ) 

OPEN (UNIT=12 , FILE=> best. gen') 

OPEN (UNIT=9 , FILE= , inp.genO 
C 

c nla=16 

c N=NLA/4 

n=15 


Maximum number of generations 

read (9,*) (aas(ij) , i j — 1 ,8) 
read (9 ,*) (bbs(ij ) , i j = 1 ,8) 
read(9,*) (hhl(ij) ,ij=l,8) 
read (9 ,*) (tthkl(ij) ,ij=l ,8) 


END 
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List 6: Modifications to of subroutine “panel. f’ 

subroutine panel (io , is , critlb ,ainipen ,nn,nd, aas ,bbs ,hhl , 
& tthkl) 


include "opt.inc" 
include "panel. inc" 
dimension is(nn,nd) 
c 

c ---------------- 

c UNIT 5 FOR THE INPUT DATA FILE 

c ---------------- 

open (5 ,f ile- 'pan . inp ' ) 
rewind 5 
c 

c ----------------- 

c UNIT 6 IS FOR THE OUTPUT DATA FILE 

open (6 ,f ile= 'pan . out ' ) 


clen = aas(is(io, 1)) 
cwid = bbs(is (io , 1) ) 
hi = hhl (is (io , 2) ) 
h2 = hi 
h3 = hi 

thick « tthkl(is(io,3)) 
igeo = 5 

call geom(sthl , sth2 , sth3 , igeo , thick) 

write (6 , *) J 

c write (6 ,*) 'Laminate =',lami 

write(6,*) 'Stiff ener thickness = ', thick 
write (6,*) 'Stiffener height = ' ,hl 
c 

c 


! modify 
! modify 

! modify 
! modify 


c [1] READING ALL INPUT DATA FOR A LAMINATE 

c USING SUBROUTINE ISKIN & STSKIN 

c 

call iskin (sthk , nmsk , exsk , eysk , vsk , gsk , amat sk , nsk , 
& lnosk,msk,tsk,thesk) 

c 

c call cskin ( sthk, nmsk, exsk, eysk, vsk, gsk, amat sk, nsk, 

c & lnosk,msk,tsk,thesk,lami) 


! modify 
! modify 
! modify 
! modify 
! modify 
! modify 
! modify 

! modify 


write (6 ,51) critlb (io) 
f ormat (' critlb = 'el4.7) 

write(6,*) ' 

return 

end 


51 
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List 7: Input file for optimization (inp.gen) 


6.667.5.7 1428 .5.00. 4. 444 . 4 . 444 .4.0.3. 6363 . 3 . 33333 ! a 

11.2.10.1818.8.61538.8.0. 7.4666.7.0.6.2222.5.894 ! b 

0 . 49375 , 0 . 50000 , 0 . 50625,0.51250,0. 51875 , 0 . 52500 ,0.53125,0. 53750 ! h 
0.060,0.066,0.072,0.078,0.084,0.090,0.096,0.102 ! t 


List 8: Input file for problem parameters (pan.inp) 


0.036 ! skin thickness 
1 ! Number of material 

I,20.2e6,1.9e6,0.3,0.73e6 ! material No., Ell, E22, vl2, G12 
6 ! number of plies 

1.1.0. 006.60.0 ! layer No., material no., ply thickness, theta 

2.1.0. 006.0.0 

3.1.0. 006,-60.0 

4.1.0. 006,-60.0 

5.1.0. 006.0.0 

6.1.0. 006.60.0 


1 

1,20. 2e6 , 1 . 9e6 , 0 
1 

1 , 1 , 0.0 
1 

1,20. 2e6 , 1 . 9e6 ,0 
1 

1 , 1 , 0.0 
1 

1,20. 2e6 , 1 . 9e6 , 0 
1 

1 , 1 , 0.0 

0 . 0 , 0 . 0 , 1 . 0el2 

20.0. 0.0.1. 0el2 

20.0. 56.0.1. 0el2 
0.0,56.0,1. 0el2 

1 . 0 . 1.0 
0 , 1 , 0,1 
1 , 1 , 1,1 
0 , 0 , 0,0 
0 , 0 , 0,0 
55 

400.0,0.0,0.0 


3,0 .73e6 


3,0 .73e6 


3,0 .73e6 


Number of material (X-stiff ener) 
material No., Ell, E22, vl2, G12 
number of plies 
layer No., material no., theta 
Number of material (Y-stiff ener) 
material No., Ell, E22, vl2, G12 
number of plies 
layer No., material no., theta 
Number of material (D-stiff ener) 
material No., Ell, E22, vl2, G12 
number of plies 
layer No., material no., theta 
1 
2 
3 


x,y,r of node 
x,y,r of node 
x,y,r of node 
x,y,r of node 
B.c's of U 
B . c ’ s of V 
B.c’s of W 
B.c ; s of Px 
B.c’s of Py 
# of modes considered 
Nx, Ny, Nxy 
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List 9: Output file containing optimal designs (best.gen) 


POPULATION SIZE= 8 CROSSOVER PROB. =1.000 
MUTATION PROB .=0.100 PERMUTATION PROB . =0 . 950 
BESTS DESIGNS AFTER 211 EVALUATIONS OF THE OF 
FS= 0 . 7669355E+03 GENERATION 3 15 
8 3 1 

FS= 0 . 1063876E+03 GENERATION 3 1 
8 1 1 

FS= 0 . 6040214E+01 GENERATION 3 1 
8 5 5 

FS= 0 . 8574801E+00 GENERATION 3 1 
7 8 5 

FS= 0 . 5234697E+00 GENERATION 3 1 
6 4 4 

FS= 0. 1314688E+00 GENERATION 3 1 
4 8 3 

FS= 0 . 4932792E-01 GENERATION 3 1 
2 6 1 

FS= 0 . 3394876E-01 GENERATION 3 1 
2 3 8 

FS= 0 . 3239560E-01 GENERATION 3 1 
113 

CPU TIME 3 516.6477 
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List 10: Output file containing convergence history (on. gen) 

population size= 8 crossover prob. =1.000 
mutation prob. =0.100 permutation prob. =0.950 
Iteration Average Best 

1 0 .4257929E+02 0 . 1063876E+03 

2 0 . 5292514E+02 0 . 1063876E+03 

3 0 . 8050877E+02 0 . 1063876E+03 

4 0 . 6810519E+02 0 . 1063876E+03 

5 0 . 1045733E+03 0 . 1063876E+03 

6 0 . 9221518E+02 0 . 1063876E+03 

7 0 . 9460207E+02 0 . 1063876E+03 

8 0 . 8101215E+02 0 . 1063876E+03 

9 0 . 7896147E+02 0 . 1063876E+03 

10 0 . 6563631E+02 0 . 1063876E+03 

11 0 . 5746792E+02 0 . 1063876E+03 

12 0 . 8395988E+02 0 . 1063876E+03 

13 0 . 6175915E+02 0 . 1063876E+03 

14 0 . 5932686E+02 0 . 1063876E+03 

15 0. 1258099E+03 0 . 7669355E+03 

16 0 . 2046968E+03 0 . 7669355E+03 

17 0 . 2256355E+03 0 . 7669355E+03 

18 0 . 1188792E+03 0 . 7669355E+03 

19 0 . 2204189E+03 0 . 7669355E+03 

20 0 . 2373301E+03 0 . 7669355E+03 

21 0.4233659E+03 0 . 7669355E+03 

22 0 . 2582285E+03 0 . 7669355E+03 

23 0 . 6843670E+03 0 . 7669355E+03 

24 0 . 3967787E+03 0 . 7669355E+03 

25 0 . 3142109E+03 0 . 7669355E+03 

26 0 . 1357727E+03 0 . 7669355E+03 

27 0 . 1403906E+03 0 . 7669355E+03 

28 0 . 2964355E+03 0 . 7669355E+03 

29 0.3865364E+03 0 . 7669355E+03 

FINAL POPULATION AFTER 211 EVALUATIONS OF THE O.F. 
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List 11: Information about first optimal design from “pan. out” 


Stiffener thickness = 5 . 9999999999999998E-02 

Stiffener height = 0.5062500000000000 

SKIN LAMINATE 

LAYER MATERIAL THK ORIENTATION 

No. No. (in) (deg) 

1 1 0.0600 0.0000 
STIFFENER ORIENTATION (deg) = 29.49017008775938 

UNIT CELL LENGTH (in) = 3.333330000000000 

UNIT CELL WIDTH (in) * 5.894000000000000 

« ONLY AXIAL & DIAGONAL STIFFENERS » 

Stiffening Parameter (X) = 0.7274678814806318 

Stiffening Parameter (Y) = 0 . OOOOOOOOOOOOOOOOE+OO 

Stiffening Parameter (D) = 0.1996764909471687 

X-stif f ener 

znn = -0.1317762672296654 

zstar = -1 . 3303185398854778E-02 

ystar = 1.616720000000000 

nstep = 56 

D-stiff ener 

znn = -0.1312237003831008 

zstar a -1 . 5189387389525246E-02 
ystar = 3.292386964576168 

nstep = 99 

Global Lamda = 1.034195740432642 

(Nx)_sk = 238.2547676115182 

(Ny)_sk = 0.0000000000000000E+00 

(Nxy)_sk = 0.0000000000000000E+00 

(Nx)_l = 912.3437009050962 

Ndx.l = 108.8378075805412 

Ndx_2 = 0. OOOOOOOOOOOOOOOOE+OO 

(Ndxy) _x = 0.0000000000000000E+00 

(Ndxy)_y = 0.0000000000000000E+00 

skilam = 1.002140728935756 

riblx = 2.002348947521088 

ribld = 33.41810771983094 

Volume (lbs/ft"2) (rho=0 . 0570)= 0.5519452823098701 

critlb = 0 . 7669355E+03 
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List 12: Information about second optimal design from “pan. out” 


Stiffener thickness = 5 . 9999999999999998E-02 

Stiffener height = 0.4937500000000000 

SKIN LAMINATE 

LAYER MATERIAL THK ORIENTATION 

No. No. (in) (deg) 

1 1 0.0600 0.0000 
STIFFENER ORIENTATION (deg) = 29.49017008775938 

UNIT CELL LENGTH (in) = 3.333330000000000 

UNIT CELL WIDTH (in) = 5.894000000000000 

« ONLY AXIAL k DIAGONAL STIFFENERS » 

Stiffening Parameter (X) = 0.7095057115675297 

Stiffening Parameter (Y) = 0 . 0000000000000000E+00 

Stiffening Parameter (D) - 0.1947462072200780 

X-stiff ener 

znn = -0.1278212016698681 

zstar = -1 . 3189737836 119143E-02 

ystar = 1.587850000000000 

nstep = 55 

D-stiff ener 

znn = -0.1272756072256236 

zstar = -1.4851388047712849E-02 

ystar = 3.292386964576168 

nstep = 99 

Global Lamda = 0.9691460449089669 

(Nx)_sk = 225.6313101069534 

(Ny)_sk = 0.0000000000000000E+00 

(Nxy)_sk = O.OOOOOOOOOOOOOOOOE+OO 
(Nx) _1 = 864.0049748708219 

Ndx_I = 103.0712516679092 

Ndx_2 = 0.0000000000000000E+00 

(Ndxy) _x = O.OOOOOOOOOOOOOOOOE+OO 
(Ndxy) _y = O.OOOOOOOOOOOOOOOOE+OO 
skilam = 1.058207774362076 

riblx = 2.147540224369275 

ribld = 35.85353200238720 

Volume (lbs/ft~2) (rho=0 . 0570)= 0.5456130037343178 

critlb = 0 . 1063876E+03 
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Chapter 1 
INTRODUCTION 


The theoretical foundations and user instructions for a FORTRAN code for the 
buckling analysis for anisotropic variable curvature panels are presented. The variable 
curvature panel is assumed to consists of two or more panels of constant curvature 
where each panel may have a different curvature. Bezier polynomials are used as 
Ritz functions. Displacement (C°), and slope (C 1 ) continuities between segments are 
imposed by manipulation of the Bezier control points. A first-order shear- deformation 
theory is used in the buckling formulation. 

Chapter 2 gives an account of the theoretical foundations of the buckling analysis. 
Instructions for setting up input files for the FORTRAN code are given in Chapter 3 
for three structure cases. Examples of input files for the structural cases considered 
are also given. Finite element results of the structural cases considered are given so 
as to provide a comparison between results from the present analysis and those of 
finite element simulations. 
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Chapter 2 
THEORETICAL 
FOUNDATIONS 


The present analysis method for buckling of anisotropic shells with variable cur- 
vature uses a segment approach where displacement fields within each segment are 
represented by Bezier polynomials and a first-order, shear-deformation theory is used. 
In general, segments can be used in both axial and circumferential directions, however 
the present implementation considers only segments in the circumferential direction. 
This restriction is based on the typical frame spacing for fuselage structures and on 
limiting general buckling to frames rather than across frames. Continuity of displace- 
ment at the junctures of adjacent segments are imposed using C° and C 1 conditions 
obtained from the properties of the Bezier control points ([1]). The shell with variable 
curvature is assumed to consist of two or more curved panels of constant curvature 
which is representative of fuselage or wing structures. 

The following sections 

• Geometry of variable curvature panel. 

• Bezier polynomials. 

• Continuities along segment junctures. 

• Minimum potential energy 


describe the formulation of the buckling analysis. 
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2.1 GEOMETRY OF VARIABLE CURVATURE PANEL 

The coordinate system and the displacement directions for a noncircular shell 
is shown in Figure 1. Any point in the wall of the shell is specified by means of 
curvilinear coordinate system x , y and z, where x is the axial coordinate fixed to 
mid-surface, y is the circumferential coordinate which follows the median line of the 
transverse cross section, and z is the radial coordinate normal to both x and y. The 
noncircular shell is assumed to consist of two or more segments in the circumferential 
direction each of constant radius. The normal and tangent vectors of the two segments 
at a juncture are equal as shown in Figure 1, where H\ = n 2 and ti = t 2 . 


2.2 BEZIER POLYNOMIALS 


Bezier polynomials are used in the axial and circumferential directions to rep- 
resent the displacement fields. The Bezier polynomial in terms of an independent 
variable is given by 




nl 


i - 1 


(i — 1)! (n — i 4- 1)! 


(u-l) 


n—i+1 


(i) 


where n denotes the order of the polynomial and 0 < v < 1. For a Bezier polynomial 
of order n, there are (n + 1) control points. The values of the control points will 
determine the variation of /,(rc, u) within the interval of 0 < v < 1. Any point on the 
surface of the segment is given by a parametric function in two variables of the form 

Prsi^l) = ( 2 ) 

r = 1 


where the coordinates £ and tj are defined as 


£ = x / L 

v = (y ~ Vi) / (2/*+i - yi) ( 3 ) 


with 0 < £, 7 ] < 1, X and Y are the number of control points in the axial and circum- 
ferential direction respectively, and q rs are the Bezier control points or coefficients. 
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The displacement vector can be written as 


' Uo ' 


' Prs 

0 

0 

0 

0 - 


qirs 

Vo 


0 

Prs 

0 

0 

0 


q2rs 

w 

> = 

0 

0 

Prs 

0 

0 

< 

q3rs 

fix 


0 

0 

0 

p 

1 rs 

0 


q4rs 

- fiy ' 

i 

. 0 

0 

0 

0 

Prs - 

3 

- q5rs 


where uq and Vq are the axial and transverse membrane displacements, respectively, 
w is the normal displacement, and (f> x and (j) y are the bending rotations. There are 
five degrees of freedom (N DOF=5) per control point and the range of subscript j 
is 1, 2, 3, ... (XY). The control points for each degree of freedom can be used to 
impose boundary conditions on each degree of freedom on each segment. 


2.3 CONTINUITIES ALONG SEGMENT JUNCTURES 

Continuity of displacement functions along segment junctures are obtained by 
using the relations between control points of the adjacent segment based on C° and 
C l continuities. Figure 2 shows two adjacent segments and the control points that 
are involved in the C° and C 1 continuities for the case of eleven control points in the 
axial direction and six control points in the transverse directions, i.e., X — 11 and V 
= 6. The control points shown in Figure 2 are for one degree of freedom and therefore 
subscript 1,2, ...,5 in q \ rs , q 2 T3 , • ••, Qsrs of Equation (4) have been dropped. In the I th 
segment, control points qke and qk 5 are related to control points qk\ and qk 2 of the 
(I + l) th segment, where k = 1,2, ..., 11 according to 


qi be 
qke 


= qk 1 for C° continuity 

= 5j qk l + pH for C 1 continuity 

Si + Si + 1 


(5) 


where 5/ and 67+1 are the width of the I th and (/ + l) th segment, respectively. Using 
these conditions the unknowns qk 1 and qk 2 are expressed in terms of qk 5 and qke,, which 
the master control points. 


The procedure for slaving qk 1 and qk 2 to qk 5 and qke, is demonstrated below 
considering only control points with subscript k and only one degree of freedom. 
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Equation (5) can be written as 


C +l) - si? 
«r , 4 / 2 +1) + « - si? 


( 6 ) 


where superscript (/ + 1) and (I) have been added to denote segment (/+ 1) and (/), 
respectively. Equation (6) can be written in matrix form as 

s£l +1) 
sll +1) 

si? 
si? 


0 

7H 

l k2 


0 a\ ,+1) 


a 


(/) 

k5 


-1 

-1 


sir 


si? 


si? 


> = 


( 7 ) 


In matrix notation form Equation (7) can be written as 




f {De} | 

’ [Cel [Cr] 

to] ; 

1 {Dr} J 


f {De} 1 _ 

[C e r] 2X 2 ^2x8 

l {Dr} J 

ffllOxiO 


- [ 0 ] ( 8 ) 

Using Equation (8), a transformation matrix [T] is developed, and 

{ D r } = [T] { D r } (9) 

where [C er ] = — [C e ] -1 [C r ]. These matrices have to be set up for k — 1, 2, 3, 

11, and for every degree of freedom. The complete matrices [C e r] 7 [Ce]i and 
[C r ] including eacn k and degree of freedom contains some populated blocks and 
matrix multiplication and inversion can be performed by manipulating the blocks for 
computational efficiency. The modified stiffness matrix [K] is given by 

[K] = [Tf[K][T] (10) 

The above matrix multiplication is performed by manipulating the populated blocks 
in matrix [T] for computational efficiency. 
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Since the buckling analysis involves first-order, shear-deformation, only C° con- 
tinuity is required in the variational formulation. However the advantage of also 
imposing C 1 continuity is not only to obtain a more accurate analysis but also to 
reduce the size of the stiffness and geometric stiffness matrices when a larger number 
of segments are used to represent the shell that is being analyzed. Table 1 shows 
the size of the matrices as indicated by the parameter ISIZE with the number of 
segments for different conditions of continuities when X = 11 and Y — 6. ISIZE is 
the number of terms in the matrix and also the number of unknown coefficients to 
be determined. If the segments are joined to approximate a closed shell, the size of 
the final matrices is less than that for a panel. The size of the stiffness and geometric 
stiffness matrices after assembly is given by the expression 


ISIZE = NDOF xX xY x NSEG - NDOF x2xX x MJOIN (11) 

for C° and C 1 continuities, where ISIZE is the matrix size or number of unknown 
coefficients, NSEG is the number of segments, and MJOIN = NSEG— 1 for a panel 
and MJOIN — NSEG for a closed shell, i.e., MJOIN is the number of junctures. 


2.4 MINIMUM POTENTIAL ENERGY PRINCIPLE 

The linear stiffness matrices are derived from the strain energy which is given by 


U 


r _ 

Aij 

B^ 

0 ' 

L {e} 



0 


_ 0 

0 

►Q 

1 


{e}dA 


( 12 ) 


where A{j is the extensional stiffness coefficient matrix, Bij is the coupling stiffness 
coefficient matrix, Dij is the bending stiffness coefficient matrix and C vq is the trans- 
verse shear stiffness coefficient matrix. The strain vector is {e} and given by 


{e} — {e x 6 y 7 X y K x Ky Kxy Ixz lyz} 


( 13 ) 


The strain-displacement relations are 

du 0 


e? = 


4 = 


dx 

dv 0 wp 

dy R 
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/C-i, — 




7 ° 

Ixz 


lyz 


du 0 dv 0 
dy + dx 
d(f) x 
dx 

d(j) y 
dy 
d(j> x 
dy 

(f>x 4" 


d (j^ C 2 _ ( dvo _ du 0 

dx 2 R dx dy 

dw 0 
dx 


, .dwo r vp 
^ + ^~ C 'R 


( 14 ) 


Here C\ and C 2 are “tracer” coefficients used to implement different strain-displace 
-ment relations or shell theories. Accordingly when C\ = C 2 — 1, the first approxima- 
tion of Sanders-Koiter shell theory [2, 3] is obtained, and when C\ = 1, C 2 = 0, Love’s 
shell theory [4] including transverse shear deformations is obtained. Finally, when C\ 
— 0 and C 2 = 0, Donnell’s shell theory [5] including transverse shear deformation is 
obtained. 


The geometric stiffness matrix is derived from the work done, (Wj), by the 
applied prebuckling loading and is given by 


W d 


-h 


Nx^xNL 4“ N y€yN L 4“ ^xylxyNL )dA 


(15) 


where the nonlinear strain components are 

(7 Xy^NL VOiy (vpyy 4“ — ) VQ, x Uq jX 4 “ V), X (W,y ) (fd) 


In the present analysis, the applied prebuckling loading is prescribed as a uniform in- 
plane stress state. The linear stiffness and geometric stiffness matrices are developed 
using analytical integration rather than numerical integration for computational effi- 
ciency. Finally, an eigenvalue problem is solved for determining the critical buckling 
load. 



Chapter 3 

USER INSTRUCTIONS 


User instructions are provided for the FORTRAN code for the buckling analysis 
of variable curvature shell. The source code is found in directory “variable”. A make 
file is used to link all the subroutines. The input file to the code is “bez.inp” and 
the out put file is “b.out”. The executable for the code is “run” as specified in the 
makefile. 

Results are presented for a composite cylindrical panel subjected to axial com- 
pression, a composite wing leading-edge panel subjected to combined axial compres- 
sion and shear, and a composite and isotropic non-circular fuselage. Sanders-Koiter 
([2, 3]) shell theory is used in these examples. Buckling loads from the present analy- 
sis are compared with those obtained from the STAGS ([6]) finite element code. The 
STAGS finite element model consists of the 410 element, curved surfaces are approxi- 
mated as an assembly of flat surfaces, and the formulation of the 410 element is based 
on the classical laminated plate theory. The nominal ply mechanical properties for 
the composite material used are: Eu — 13.75 Msi; E 2 2 = 1.03 Msi; G\ 2 =G\ 3 =G 2 3 = 
0.420 Msi and Vi 2 = 0.250, The laminate ply stacking sequence is [±45/0/90/ ± 45 ] s 
with equal ply thicknesses for each of the different laminate thicknesses. For the 
isotropic material, the mechanical properties are Eu — E 22 = 10.0 Msi; v\ 2 = 0.3, 
and Gi 2 ~Gi 3 =G 2 3 — Eu / 2(1 ± ^ 12 ). 

3.1 APPLICATION OF BOUNDARY CONDITIONS 

Boundary conditions can be imposed along £ = 0,1, r/ = 0,1 (Figure 2) and also 
at any control point by the subroutine “ibcs.f”. For example, if w=0 at £ = 1 and w 
— 0 at 7/ = 0 on segment 1, then these constraints are entered as 


9 


2 ! Number of line constraints 

1,1, 1,1 ! coordinate No., coordinate value, segment No., Dof No. 

2,0, 1,3 

Characters after the “!” are only comments. Coordinate No. can be either 1 or 2, 
coordinate No. = 1 for x and coordinate No. = 2 for y. Coordinate value = 0 or 1, 
coordinate value = 0 for £ or rj = 0, and coordinate value = 1 for £ or y = 1. Segment 
No. is the segment No. where the constraint is applied. Dof No. can be 1, 2, 3, 4, or 
5 and u= 1, v=2, w=3, (f) x ~ 4, and <f> y = 5. 

Point constraints can be applied by specifying the control point associated with 
the degree of freedom of concern. For example if control point q rs , where r = 1,2, 
3, ... ,11, and 5=1,2, 3, ... , 6, (see Figure 2) of segment ( I -f 1) is to have w= 0 
constraint, then the number of point constraint is one and the entry format to be 
entered for q rs with w = 0 is 

No. of point constraints 

segment No., control point No., Dof No 

The segment No. is (7+1), the control point No. for q rs is [(r — 1) x 6 + s], and Dof 
No. is as described above. If there is no point constraint then the number of point 
constraint is entered as zero. 

3.2 CYLINDRICAL PANEL 

The first structure analyzed is a semi-circular laminated composite (a = 180°) 
cylindrical panel 22.0-in. long, and with a radius of 40.0 inches as shown in Figure 3. 
The simply-support boundary conditions are also shown in Figure 3. The cylindrical 
panel is modeled with five curved segments in the present analysis while the STAGS 
finite element modeled consists of a mesh of 20 x40 elements in the axial and transverse 
direction, respectively. Table 2 shows the results for the curved panel subjected to 
axial compression load for different thicknesses, List 1 shows the input data for the 
problem where characters after the “!” are only comments and need not be included 
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in the input file. The input file is “variable/bcpan5.inp” on the compact disk and the 
file has to be named “bez.inp” for execution. 

The results in Table 2 suggest that for t = 0.072 in. the present analysis result is 

4.3 percent greater than the STAGS result, while for t = 0.144 in. and 0.216 in. the 
STAGS analysis results is 1.4 percent above that of the present analysis. The STAGS 
results are above that of the present analysis for t — 0.144 in. and 0.216 in. since 
these two panels are thicker and hence transverse shear deformation effects have an 
influence on buckling load. 

3.3 WING LEADING-EDGE PANEL 

The wing leading-edge panel is shown in Figure 4. It consists of three curved 
segments of radii 50.0 in., 6.136 in., and 50.0 in., respectively, and is 26.0-inches long, 
and has a maximum width of 32.0 inches. The boundary conditions of the panel is 
shown in Figure 4 and correspond to classical simply support conditions. Each ply 
of the laminate is 0.006-in. thick. Using the present analysis, the wing leading-edge 
panel is modeled as a combination of two segments for the 50.0-in. -radius section 
and one segment for the 6. 136-in. -radius section. The wing leading-edge panel is also 
modeled as a quarter symmetric model in the present analysis as shown in Figure 

5 with symmetry boundary conditions on two edges and classical simply support 
conditions on the other two edges. A combination of two segments for the 50.0-in.- 
radius section and one segment for the 6. 136-in. -radius section is used for the quarter 
symmetry model. The STAGS finite element model consists of the 410 shell element 
and 30x30 elements in each curved segment. 

Table 3 shows the results obtained from the present analysis for the full model 
and the quarter symmetric model for some selected combined load cases. Results from 
STAGS for the same selected combined load cases are also shown in Table 3. Figure 

6 shows the buckling load interaction curve between axial compression and positive 
shear loading. The results from the present analyses for the full model are about 4.0 
to 5.0 percent above those of STAGS except for the case of negative shear loading 
where the result from the present analysis is 7.0 percent above that of STAGS. 


11 


Better agreement with STAGS is obtained from the present method by using a 
quarter symmetric model except for the case where the shear loading is dominant as 
is seen in Figure 6. For example, there is only a very small difference between the 
results of the present analysis for the full and quarter symmetric model for the pure 
shear case. For the case where the loading condition is N x =N xy = 1, the result for 
the quarter symmetric model is about 3 percent above that of STAGS, and for the 
case of compression only, the result for the quarter symmetric model is in very good 
agreement with that of STAGS. 

The quarter symmetric model provides results in better agreement compared to 
those of STAGS for cases where the buckling mode shape is symmetric or close to 
being symmetric since fewer control points can be used to provide a more accurate 
model of the structure. The symmetric model is not recommended for loading cases 
where the shear load is predominant. The input files for this problem is given in List 
2 and 3 for the full model and quarter symmetric model, respectively. Comments 
to the input data are provided after the “!” in List 2 and 3. The input file is “vari- 
able/bw5.inp” and “variable/bwss3.inp” for the full model and quarter symmetric 
model , respectively, on the compact disk and any file has to be named “bez.inp” for 
execution. 

3.4 NON-CIRCULAR FUSELAGE 

A non-circular fuselage section is shown in Figure 7 and is 24-inches long. It 
consists of 30.0-in. -radius curved segments and flat segments, and has simply support 
boundary conditions. In both the present analysis and the STAGS analysis, a quarter 
symmetric model of the non-circular fuselage section is considered, and the boundary 
conditions are shown in Figure 8. In the present analysis, the quarter symmetric 
model of the fuselage consists of a combination of four segments, one for each flat 
and curved segment. The STAGS model consists of 410 shell element and a mesh of 
40x40 elements in each segment. 

The results are shown in Table 4 for the isotropic and laminated non-circular 
fuselage with different wall thicknesses subjected to axial compression. The buckling 
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loads for the isotropic non-circular fuselage obtained from the present analysis are less 
than 2.5 percent above those of STAGS for the different wall thicknesses considered. 
For the laminated non-circular fuselage, the buckling loads are in very good agreement 
with the those of STAGS for the wall thicknesses considered. The input file for this 
problem is given in List 4 for the isotropic case. Comments to the input data are 
provided after the “!” in List 4 The input file is “variable/bfs4.inp” on the compact 
disk and the file has to be named “bez.inp” for execution. 
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LIST OF FILES 


List 1: Input file for cylindrical panel. 


0.216 ! skin thk 

1 ! # of material 

1 . 0 . 1375e+8 , 1 . 030e6 ,0 . 25 , 0 .42e6 ! mat #, el,e2,nu,gl2 

12 ! # of ply 

1.1.0. 018.45.0 ! layer #, material #, layer thk, theta 

2.1.0. 018,-45.0 

3.1.0. 018.0.0 

4.1.0. 018.90.0 

5.1.0. 018.45.0 

6.1.0. 018,-45.0 

7.1.0. 018,-45.0 

8.1.0. 018.45.0 

9.1.0. 018.90.0 

10.1.0. 018.0.0 

11.1.0. 018,-45.0 

12.1.0. 018, 45.0 

5 ! # of segment 

22.0. 25.132741.40.0 ! length, width, radius of segment i 

22.0,25.132741,40.0 

22.0,25.132741,40.0 

22.0,25.132741,40.0 

22.0. 25.132741.40.0 

24 ! # of line constraints 

1.0. 1.2 ! Coord. #, Coord Val. , segment #, DoF # 

1.0. 1.3 

1 . 0 . 2. 2 

1 . 0 . 2. 3 

1.0. 3. 2 

1 . 0 . 3. 3 

1 >0,4,2 

1 . 0 . 4. 3 

1 . 0 . 5 . 2 

1.0. 5. 3 x=0, v=w=0 

1,1, 1,2 

1.1. 1.3 

1, 1,2,2 

1 . 1 . 2. 3 

1. 1.3.2 

1.1. 3. 3 

1,1, 4, 2 
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1 , 1 , 4, 3 

1.1. 5. 2 

1. 1.5.3 x=l, v=w=0 

2 , 0 , 1,1 

2.0. 1.3 y=0, u=w=0 

2. 1.5.1 

2. 1.5. 3 y=l , u=w=0 

0 ! # of point constraints 

4 ! # of joints 

2.1 ! segment 2 join to segment 1 

3.2 ! segment 3 join to segment 2 

4.3 ! segment 4 join to segment 3 

5.4 ! segment 5 join to segment 4 

1.0. 0.0. 0.0 ! Nx, Ny, Nxy 
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List 2: Input file for wing leading-edge panel (full model). 


0.072 ! skin thk 

1 ! # of material 

1 . 0 . 1375e+8 , 1 . 030e6 , 0 . 25 , 0 .42e6 ! mat #, el,e2,nu,gl2 

12 ! # of ply 

1.1.0. 006.45.0 ! layer #, material #, layer thk, theta 

2.1.0. 006,-45.0 

3.1.0. 006.0.0 

4.1.0. 006.90.0 

5.1.0. 006.45.0 

6.1.0. 006,-45.0 

7.1.0. 006,-45.0 

8.1.0. 006.45.0 

9.1.0. 006.90.0 

10.1.0. 006.0.0 

11.1.0. 006,-45.0 

12.1.0. 006, 45.0 

5 ! # of segment 

26.0. 12.606817.50.0 ! length, width, radius of segment i 

26.0. 12.606817.50.0 

26.0. 9. 037922 . 6 . 1362094 

26.0. 12.606817.50.0 

26.0. 12.606817.50.0 

24 ! # of line constraints 

1.0. 1.3 ! Coord. #, Coord Val., segment #, DoF # 

1.1. 1.3 

2. 0. 1. 3 seg. 1, w=0 

1.1. 2. 3 

1.0. 2. 3 seg* 2, w=0 

1. 1.3.3 

1.0. 3. 3 seg. 3, w=0 

1. 1.4.3 

1.0. 4. 3 seg. 4, w=0 

1.1. 5. 3 

1.0. 5. 3 

2. 1.5.3 seg. 5, w=0 

2 . 0 . 1.1 

1,0, 1,2 

1.1. 1.2 seg. 1, in-plane b.c 

1. 1.2.2 

1.0. 2. 2 seg. 2, in-plane b.c 

1. 1.3.2 

1.0. 3. 2 seg. 3, in-plane b.c 

1.1. 4. 2 

1.0. 4. 2 seg. 4, in-plane b.c 
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1. 1.5.2 

1 , 0 , 5, 2 

2. 1.5.1 seg. 5, in-plane b.c 

0 ! # of point constraints 

4 ! # of joints 

2.1 

3.2 

4.3 

5.4 

0 . 0 , 0 . 0 , 1.0 


! Nx, Ny, Nxy 
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List 3: Input file for wing leading-edge panel (quarter symmetric model). 


0.072 ! skin thk 

1 ! # of material 

1.0. 1375e+8 , 1 . 030e6 , 0.25,0. 42e6 ! mat #, el,e2,nu,gl2 

12 ! # of ply 

1.1.0. 006.45.0 ! layer #, material #, layer thk, theta 

2.1.0. 006,-45.0 

3.1.0. 006.0.0 

4.1.0. 006.90.0 

5.1.0. 006.45.0 

6.1.0. 006,-45.0 

7.1.0. 006,-45.0 

8.1.0. 006.45.0 

9.1.0. 006.90.0 

10.1.0. 006.0.0 

11.1.0. 006,-45.0 

12.1.0. 006, 45.0 

3 ! # of segment 

13.0. 12.606817.50.0 ! length, width, radius of segment i 

13.0. 12.606817.50.0 

13.0. 4.518961.6.1362094 

16 ! # of line constraints 

1.0. 1.1 ! Coord. #, Coord Val., segment #, DoF # 

1 . 0 . 2.1 

1.0. 3. 1 sym 

1.0. 1.4 

1 . 0 . 2. 4 

1.0. 3. 4 sym 

1.1. 1.3 

1.1. 2. 3 

1.1. 3. 3 w=0 

1 , 1 , 1,2 

1, 1,2,2 

1.1. 3. 2 v=0 

2.0. 1.3 

2. 0. 1.1 segment 1, u=w=0 

2. 1.3. 2 

2, 1 ,3,5 sym 

0 ! # of point constraints 

2 ! # of joints 

2,1 

3.2 

0 . 8 , 0 . 0 , 1.0 


! Nx, Ny, Nxy 
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List 4: Input file for non-circular fuselage (isotropic) 


0 . 100 ! skin thk 

1 ! # of material 

1.0. 106+8,0.106+8,0.30,3846153.8 ! mat #, el , e2 ,nu,gl2 -isotropic(Al) 

2 ! # of ply 

1.1.0. 050.0.0 ! layer #, material #, layer thk, theta 

1.1.0. 050.0.0 

4 ! # of segment 

12 . 0 . 21 . 213203 . 1 . 0e8 ! length, width, radius of segment i 

12.0. 70.685835.30.0 

12.0. 30.0.1. 0e8 

12.0. 23.561945.30.0 

16 ! # of line constraints 

1.0. 1.3 seg. 1, w=0 ! Coord. #, Coord Val., segment #, DoF # 

1.0. 2. 3 seg. 2, w=0 

1.0. 3. 3 seg* 3, w=0 

1.0. 4. 3 seg. 4, w=0 at x=0 or xi=0 

1.1. 1.4 seg. 1 

1.1. 2. 4 seg. 2 

1.1. 3. 4 seg. 3 

1.1. 4. 4 seg. 4 phi_x =0 at x=12 or xi=l 

1.1. 1.1 seg. 1, u=0 

1.1. 2.1 seg. 2, u=0 

1.1. 3.1 seg. 3, u=0 

1.1. 4.1 seg. 4, u-0 at x=12 or xi=l 

2.0. 1.2 seg. 1, v=0 at y=0 or eta=0 

2. 0. 1. 5 seg. 1, phi_y=0 at y=0 or eta=0 

2. 1.4. 2 seg. 5, v-0 at y=23, 561945 or eta=l 

2. 1.4. 5 seg. 5, phi_y=0 at y=23, 561945 or eta=l 

0 ! # of point constraints 

3 ! # of joints 

2,1 

3.2 

4.3 

1.0. 0.0. 0.0 ! Nx, Ny, Nxy 
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Table 1 Size of stiffness matrices for panel and shell for increasing number of 

segments. 
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NSEG 

ISIZE 

(Panel) 

[C°] [C°, C 1 ] 

ISIZE 
(shell) 
[C°, C 1 ] 

1 

330 

330 


2 

605 

550 

440 

3 

880 

770 

660 

4 

1155 

990 

880 

5 

1430 

1210 

1100 

6 

1705 

1430 

1320 


Table 2 Comparison of buckling loads results for composite curved panel. 


Thickness STAGS Present 
t (in.) analysis 

(lbs/in.) (lbs/in.) 


0.072 

0.144 

0.216 


374.55 390.68 

1481.08 1459.45 

3328.25 3278.86 
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Table 3 Comparison of buckling loads results for wing leading-edge panel. 


Loading 

condition 

| 

Full Model 

Quarter 

Model 

STAGS 

Present 

Analysis 

Present 

Analysis 

N x 

N X y 

(lbs/in.) 

(lbs/in.) 

(lbs/in.) 

1.0 

0.0 

301.73 

317.03 

299.93 

1.0 

0.4 

193.96 

202.57 

196.97 

1.0 

1.0 

106.34 

110.90 

108.82 

0.4 

1.0 

123.93 

129.53 

127.93 

0.0 

-1.0 

-117.13 

-125.37 

-125.60 

0.0 

1.0 

138.78 

l 

_ 

145.55 

144.64 


Table 4 Comparison of buckling loads results for non-circular fuselage. 


Thickness 

STAGS 

Present 

t (in.) 

(lbs/in.) 

Analysis 

(lbs/in.) 


Isotropic 

- 

0.080 

14.64 

14.90 

0.100 

28.18 

28.84 

0.144 

82.02 

84.13 

0.180 

157.52 

160.79 



Laminated 


0.072 

6.34 

6.32 

0.096 

14.77 

14.77 

0.120 

28.40 

28.41 

0.144 

48.39 

48.29 
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Figure 3: Geometry and boundary conditions of curved composite panel. 



32 in. 32 in. 

Rj = 50 in. R 2 = 6.136 in. 


Figure 4: Geometry, dimensions and boundary conditions of composite wing 

leading-edge panel. 
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Figure 5: Boundary conditions for quarter symmetric model of composite wing 

leading-edge panel. 



Figure 6: Buckling load interaction curve for the composite wing-leading edge panel. 
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Figure 7: Geometry of non-circular fuselage. 



Figure 8: Boundary conditions for quarter symmetric model of non-circular fuselage, 






